PL-TR-93-2070 


AD-A265  035 

iiiir . 


QUANTIFICATION  OF  mL  FOR  SMALL  EXPLOSIONS 


R.  B.  Herrmann 
K.  Hutchensen 


Department  of  Earth  and  Atmospheric  Sciences 
St.  Louis  University 
3507  Laclede  Avenue 
St.  Louis,  MO  63103 


12  March  1993 


Final  Report 

July  5, 1990  -  December  31, 1992 


Approved  for  public  release;  distribution  unlimited 


§2  5 


20  02  9 


PHILLIPS  LABORATORY 
Directorate  of  Geophysics 
AIFt^ORCE  MATERIEL  COMMAND 
HANSCOM  AIR  FORCE  BASE,  MA  01731-3010 


THIS  DOCUMENT  IS  BEST 
QUALITY  AVAILABLE.  THE  COPY 
FURNISHED  TO  DTIC  CONTAINED 
A  SIGNIFICANT  NUMBER  OF 
PAGES  WHICH  DO  NOT 
REPRODUCE  LEGIBLY. 


The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors 
and  should  not  be  interpreted  as  representing  the  official  policies,  either 
expressed  or  implied,  of  the  Air  Force  or  the  U.S.  Government. 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


JAMES  F.  LEWKOWICZ 
Contract  Manager 
Solid  Earth  Geophysics  Branch 
Earth  Sciences  Division 


1 


/ 


77 

J/(Mp6  F.  LEWKOWICZ 
Branch  Chief 

Solid  Earth  Geophysics  Branch 
Earth  Sciencs  Division 


DONALD  H.  ECKHARDT,  Director 
Earth  Sciences  Division 


This  document  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing  list, 
or  if  the  addressee  is  no  longer  employed  by  your  organization,  please  notify 
PL/TSI,  Hanscom  AFB  MA  01731-3010.  This  will  assist  us  in  maintaining  a 
current  mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on 
a  specific  document  requires  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


form  Approved 
OMB  No  0 704-0 1 88 


PuWte  rfDOrtinc  Ouraeo  »0'  I'M  Of  .ntormancn  +  e4t-ma!e<3  ro  '  *0 ur  OT<  '«OC r>e  'or  '<*v«e*v  ro  <f>V jr.;n  »*f'c*.ng  0*l*  K5wf<« 

0*th«<m0  *rx3  maintaining  the  oat  J  n«*3e<3  ano  c0mOipnng  *no  ^v.ew.ng  ;n*  collection  0»  *n»orm*t:o«  Seoa  commfnn  r*>'a'd>ng  T*<»  Ou'Oen  O'  *n?  ot**f  0»  t^'» 

collection  of  mtofmatton.  .n<iud»ng  >ugg«M>oni  tor  rwjuong  tm»  ovi^om  to  w*sn«ogton  Hp»eou*nen  Vef*icei  Directorate  for  <ntc"nat.cn  Ooe'*t  onj*no  *eo©n*,  Wt5iefte^v>r 
Highway  Suite  ’204  Arlington  VA  22J02-4J02  and  to  tne  Offi<e  O*  Management  *ng  Suoge’  **©e**wrOr».  «Wu<t*on  Project  iQ/CS-O’SBf  /va»ningtcr  D<  205C3 


1.  AGENCY  USE  ONLY  (Leave  bian *) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

12  ; lurch  1993 


Quantification  of  m,  for  Small  Explosions 

hg 


6.  AUTHOR(S) 

R.  B.  Herrmann 
K.  Hutchensen 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Earth  and  Atmospheric  Sciences 
Saint  Louis  University 
3507  Laclede  Avenue 
St.  Louis,  MO  63103 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Phillips  Laboratory 
29  Randolph  Road 
Hanscom  A IB ,  MA  01731-3010 

Contract  Manager:  James  Lewkowicz /GPEH _ 


11.  SUPPLEMENTARY  NOTES 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  public  release; 
Distribution  unlimited 


5.  FUNDING  NUMBERS 

PL  62L/ 11 

Pi:  7600  ii»  09  V.T  BC 
font rac t  F 1 962U-90-K-0G4 


8.  PERFORMING  ORGANIZATION  ■ 

REPORT  NUMBER  j 


10.  SPONSORING ,  MONITORING 
AGENCY  REPORT  NUMBER 


PL-TR-9 3-2070 


12b.  DISTRIBUTION  code 


13.  ABSTRACT  (Maximum  200  words) 

This  report  focuses  on  regional  wave  propagation  with  application  to  the  dis¬ 
crimination  problem  between  explosions  and  earthquakes.  Since  discrimination 
is  not  trivial  for  small  events,  understanding  of  seismic  wave  propagation  and 
seismic  sources  is  essential.  The  report  begins  by  completely  presenting 
seismic  wave  propagation  theory  for  plane  layers  in  an  isotropic  media. 

Topics  such  as  propagator  matrix  stability,  wavenumber  integration  and 
periodicity  of  the  discrete  Fourier  transform  are  addressed.  The  next  topic 
discusses  the  differences  between  small  explosions  and  earthquakes  of  the 
same  magnitude.  The  investigation  indicates  that  chemical  mining  explosions 
have  low  excitation  of  high  frequencies  compared  to  point  sources,  such  as 
earthquakes.  Finally,  the  use  of  the  short-period  surface  wave  as  a  depth 
indicator  is  discussed. 


14.  SU8JECT  TERMS 


Explosions,  Earthquakes,  Seismic  Wave  Theory,  Seismic  Moment 


15  NUMBER  OF  PAGES 

104 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

Unclassified  Unclassified _ Unc  1  assi  f  ied _  SAR 


NSN  7540-01-280-5500  Standard  form  298  (Rev  2-89! 

Pr«<MOCd  bv  AivSi  5?a  £39  1 8 


I 


TABLE  OF  CONTENTS 


Preface .  v 

Elastic  Wave  Green’s  Functions  for  Isotropic  Layered  Media .  1 

Introduction .  1 

Functional  Form  of  the  Solution .  2 

Wholespace  Solution .  7 

Halfspace  Solution .  12 

Plane  Layered  Media .  15 

Liquid  Layer .  22 

Time  Histories .  29 

Wavenumber  Integration .  34 

Numerical  issues  in  Evaluating  Propagator  Matrices .  44 

References .  46 

Appendix .  48 

Spectral  Examination  of  the  16  June  1992  Earthquake  and 

Quarry  Blast  near  Evansville,  Indiana .  57 

Introduction .  57 

Description  of  the  Study  Area .  60 

Data  Processing  and  Results .  60 

Discussion  and  Conclusions .  76 

Acknowledgments .  79 

References . 79 

Discourse  on  the  Use  of  Short  Period  Surface  Waves  as  a 

Depth  Discriminant .  85 

Introduction .  85 

Depth  Constraints . 86 

Discussion .  88 

References .  89 


Aaosssion  For  /  j 

Nii  j 

DTiC 

J-.:  •  1 

T  *  '■  j 

< une p d 
i  icaAion _ 

& 

□ 

□ 

Pv  .. 

Dintr 

Aval 

i-lbillty  Codes 

Dlst 

Avail  and/or 
Special 

Mi 

Preface 

This  report  focuses  on  regional  wave  propagation  with  application  to 
the  discrimination  problem  between  explosions  and  earthquakes.  Since  dis¬ 
crimination  is  not  trivial  for  small  events,  no  one  single  technique  will  resolve 
the  problem.  This  implies  that  one  must  understand  seismic  wave  propaga¬ 
tion  in  all  its  complexity.  There  are  limits  to  this  both  in  computation  and 
theory  and  in  knowledge  of  true  three  dimensional  real  earth  model. 

As  a  step  toward  understanding  wave  propagation,  the  first  part  of  this 
report  presents  a  complete  review  of  wave  propagation  in  layered,  isotropic, 
homogeneous  media.  Analytic  solutions  for  wave  generation  in  a  wholespace 
due  to  moment  tensor  and  point  forces  is  given  in  the  r-co  domain,  where  r  is 
the  radial  distance  between  the  source  and  receiver.  This  is  followed  by 
explicit  expressions  for  the  solution  to  the  same  problem  in  the  k-w  domain, 
where  k  is  wavenumber.  The  reason  for  these  two  representations  is  to  per¬ 
mit  validation  of  numerical  integration  techniques  used  to  obtain  time  histo¬ 
ries  as  a  function  of  distance.  The  solutions  in  the  k  -  a>  domain  are  extended 
for  the  halfspace  and  a  generalized  layered  medium.  The  solution  to  the  lay¬ 
ered  media  problem  is  given  in  terms  of  propagator  matrices.  The  general 
solution  permits  a  variety  of  boundary  conditions  at  the  top  and  bottom  of  the 
layered  stack  that  permits  the  same  methodology  to  be  used  for  reflectivity 
computations  of  P-wave  first  arrivals.  The  subjects  of  wavenumber  integra¬ 
tion,  propagator  matrix  stability  and  consequences  of  using  discrete  Fourier 
transforms  are  also  discussed.  The  algorithms  presented  and  the  tutorial  on 
their  use  are  being  used  to  understand  regional  wave  propagation. 

The  second  subject  discussed  partially  addresses  the  real  problems  of 
discrimination:  that  of  distinguishing  shallow  earthquakes,  shallow  delayed 
chemical  explosions  and  shallow  point  explosions.  The  real  problem  is  that 
data  sets  containing  all  three  of  these  phenomena  are  not  readily  available. 
This  section  compares  a  shallow  earthquake  and  a  delayed  explosion  of  simi¬ 
lar  magnitude  near  Evansville,  Indiana  in  1992.  Differences  in  the  two 
sources  are  seen  in  the  P-wave  spectra,  where  the  strip  mining  explosions 
have  a  lower  comer  frequency  and  significantly  less  high  frequency  signal 
above  10  Hz  results  for  the  Lg  phase  are  similar.  This  may  lead  to  a  discrimi¬ 
nant  between  spatially  distributed  chemical  explosions  and  point  earthquake 
or  explosion  sources. 

The  final  topic  touches  upon  depth  discriminant  information  contained 
in  the  Rg  wave.  Under  the  conditions  of  a  relatively  uniform  waveguide  and 
with  broadband  recording,  the  short  period  surface  wave  may  be  a  useful 
quantifiable  indicator  of  depth.  The  presence  or  absence  of  a  short  period 
surface  wave  does  not  give  any  direct  depth  information,  unless  additional 


independent  information  is  available,  such  as  knowledge  of  regional  velocity 
and  Q  models  and  absolute  source  size.  Given  these  constraints  a  family  of 
curves  can  be  created  to  indicate  the  moment  required  to  generate  an 
observed  signal  of  a  given  amplitude  as  a  function  of  depth.  If  a  suite  of 
curves,  generated  under  other  assumptions  of  earth  structure,  is  used,  a 
"fuzzy"  estimate  of  source  depth  can  be  made  simply  on  the  basis  of  whether 
a  surface  wave  of  a  given  amplitude  is  observed. 

The  implication  of  this  research  to  the  task  of  event  discrimination  is 
that  a  systematic  attempt  has  been  made  to  develop  the  tools  for  understand¬ 
ing  regional  wave  propagation  and  for  their  use  in  analyzing  events  that  are 
not  confidently  rejected  from  consideration  on  the  basis  of  exhibiting  strong 
earthquake  or  delayed  explosion  characteristics.  To  classify  these  remaining 
events,  the  entire  recorded  waveform  must  be  understood.  From  the  theoreti¬ 
cal  point  of  view,  broadband  signals  above  the  noise  level  are  preferred  to 
accomplish  this. 
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ELASTIC  WAVE  GREEN’S  FUNCTIONS  FOR  ISOTROPIC  LAYERED  MEDIA 

Robert  B.  Herrmann 
ABSTRACT 

The  Haskell  (1964)  formulation  for  the  elastic  wave  field  due  to 
point  sources  in  a  plane  layered  medium  is  revised  for  the  pur¬ 
pose  of  generalizing  the  boundary  conditions  at  the  top  and  bot¬ 
tom  of  the  layer  stack  and  for  permitting  the  receiver  to  be  at 
any  depth  in  the  medium.  Because  of  the  usefulness  and  also 
complexity  of  the  problem,  explicit  solutions  are  given  for 
wholespace  and  halfspace  problems.  The  numerical  implementa¬ 
tion  is  discussed  in  detail,  demonstrating  needed  tricks  to 
ensure  quality,  noise-free  solutions. 

INTRODUCTION 

In  testing  seismogram  synthesis  programs,  it  is  necessary  to  compare 
the  program  results  to  analytic  solutions  or  to  previous  numerical  solutions. 

This  has  led  to  several  papers  by  the  author  and  students  at  Saint  Louis  Uni¬ 
versity  (Wang  and  Herrmann.  1980;  Herrmann  and  Wang,  1985;  Herrmann 
and  Mandal,  1986).  These  solutions  were  for  center  of  expansion  and  double 
couple  sources  in  a  wholespace  and  in  a  halfspace  observed  at  the  free  sur¬ 
face.  Recently  the  author  found  it  necessary  to  consider  the  case  of  a  buried 
receiver,  for  which  a  stable  algorithm  was  required.  In  addition  closed  form 
solutions  for  point  forces  were  also  required  because  of  their  use  as  sources  in 
seismic  exploration  sources  and  because  of  interest  in  the  non-isotropic  signal 
associated  with  underground  explosions.  Because  of  the  author’s  familiarity 
with  the  Haskell  (1963,  1964)  papers,  this  presentation  will  continue  using 
that  notation  in  a  cylindrical  coordinate  system  because  of  the  assumption  of 
isotropic,  laterally  homogeneous  layering.  Other  representations  of  the  solu¬ 
tions  are  possible,  especially  the  reflectivity  coefficients  of  Kennett  (1983). 

Given  the  assumption  of  a  laterally  homogeneous  isotropic  medium  and 
because  the  Green’s  functions  will  be  given  for  point  sources,  it  is  convenient 
to  construct  the  solution  in  a  cylindrical  coordinate  system  (r,<p,z)  and  in 
terms  of  angular  frequency.  The  specific  definition  of  the  Fourier  transform, 

H(<»),  of  the  time  series  h(t)  is  taken  to  be 

oo 

H (a>)=  J  h(t)e"ia,tdt 

—  oo 

with  the  inverse  transform  defined  as 

h(t)  =  —  [  H(oi)ei<”td/y 
2k  J 
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FUNCTIONAL  FORM  OF  THE  SOLUTION 

Haskell  (1963)  built  the  solution  for  the  displacement  field  due  to  point 
couples  in  a  wholespace  by  starting  with  the  analytic  solution  for  the  dis¬ 
placement  field  due  to  a  point  force  given  in  a  cartesian  coordinate  system. 
Solutions  for  point  single  couples  were  obtained,  and  the  solution  was  cast 
into  a  cylindrical  coordinate  system,  through  the  use  of  partial  derivatives  of 
the  Sommerfeld  integral.  Haskell  (1964)  extended  the  Haskell  (1963)  work  to 
a  layered  halfspace,  including  double-couple,  dipole  and  point  forces.  The 
equations  below  cast  the  Haskell  (1963)  derivations  into  the  Green’s  func¬ 
tions  for  dislocation  and  explosive  sources,  originally  given  by  Herrmann  and 
Wang  (1985).  The  Green’s  functions  are  defined  as  follow: 


CM 

ZDD  =  J  Fx(k,  e>)Jo(kr)dk 

0 

(la) 

oo 

RDD  =  -  J  F2(k,  <u)J1(kr)kdk 

0 

(lb) 

oo 

ZDS  =  |  F3(k,  <w)J1(kr)dk 

0 

(lc) 

oo 

RDS  =  |  F4(k,  <y)J0(kr)kdk 

0 

(Id) 

1  °° 

--  f  [F4(k,  co)  +  F^k.fi^Jjikrldk 
r  J 

0 

oo 

TDS=  fF„(k,<y)J0(kr)kdk 

0 

(le) 

1  00 

--  f[F4(k,<w)  +  F13(k,  cy)]  J^krldk 
r  J 

0 

oo 

ZSS=  J  F6(k,<y)J2(kr)dk 

0 

(10 

oo 

RSS  =  J  Fe(k,  <w)J|(kr)kdk 

0 

dg) 
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~  -  J  [F6(k ,a>)  +  F  14(k,  <y)]  J2(kr)dk 
r  o 

w 

TSS  =  |  F14(k,  <a)Ji(kr)kdk  ( lh) 

o 

-  -  J  [F6(k,<y)  +  F  i4(k,  a))]  J2(kr)dk 
r  0 


oo 

ZEP=  jF7(k,^)Jo(kr)dk  (li) 

o 

oo 

REP  --  j  Fg(k,  <y)J1(kr)kdk  (lj) 

o 

oo 

ZVF  =  J  F9(k,  «y)Jo(kr)dk  (lk) 

o 

oo 

RVF  =  -  J  F  10(k,  ©VWkrikdk  ( 11) 

0 

oo 

ZHF  =  J  Fu(k,  <w)J1(kr)dk  (lm) 

o 

oo 

RHF  =  J  FI2(k,  <y)J0(kr)kdk  ( In) 

o 


--  J[F12(k,®)  +  F^k.^WkrWk 

o 

oo 

THF  =  J  F15(k,  )J0(kr)kdk  ( lo) 

o 

J[F12(k,^)  +  F  15(k, <y)]  J j(kr)dk 

r  o 
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(Ip) 


PEP  =  J  F16(k,  <y)J0(kr)dk 

o 

In  an  isct-jpic  medium,  an  arbitrarily  oriented  double  couple  without 
moment  so>'  ze  model  with  vector  n  =  (n1(n2,n3)  normal  to  the  fault  and  f 
=  (f1,f2)f3)  in  the  direction  of  the  dislocation  (  Haskell,  1963;  Haskell,  1964) 
has  the  following  Fourier  transformed  displacements  for  a  source  at  depth  h 
and  the  receiver  at  a  distance  r  from  the  origin  and  at  a  depth  z: 

ur(r,  z,  h,  <y)  =  ZSS  [(fxni  -  f2n2)  cos  2  0  +  (.fxn2  +  f2nx)  sin  20)  (2a) 

+  ZDS  [(fin3  +  f3ni )  cos  0  +  (f2  n3  +  f3n2)  sin  0) 

+  ZDD  [f3n3] 

ur(r,  z,h,<y)  =  RSS[(f1n1-f2n2)cos20  +  (f1n2  +  f2ni)sin20]  (2b) 

+  RDS  [(fjn3  +  f3nj )  cos  0  +  (f2n3  +  f3n2)  sin  <p\ 

+  RDD[f3n3] 

u#(r,  z,  h,  t o )  =  TSS  [(^n  x  -  f2n2)  sin  20  - (fjn2  +  f2ni)  cos  2<p]  (2c) 

+  TDS  t(fin3  +  f3n! )  sin  0  -  (f2n3  +  f3n2)  cos  0] 

The  vertical  displacement  u*  is  defined  positive  upward,  the  radial  dis¬ 
placement  is  positive  away  from  the  source,  and  the  tangential  displacement 
u#  is  positive  in  a  direction  clockwise  from  north  when  looking  downward 
from  above  the  source.  The  vectors  n  and  f  are  still  defined  in  a  local  coordi¬ 
nate  system  at  the  source  in  which  the  cartesian  x,y,z  axes  are  in  the  north 
(0  =0°),  east  (0  =90°)  and  downward  directions,  respectively.  Following  Her¬ 
rmann  (1975)  the  components  of  these  vectors  can  be  expressed  in  terms  of 
the  fault  plane  parameters  of  strike,  dip  and  rake.  The  strike,  0{,  is  mea¬ 
sured  clockwise  from  north,  the  dip,  is  measured  in  a  positive  sense  from 
the  horizontal  direction  perpendicular  to  strike,  and  the  rake,  Af,  is  measured 
on  the  fault  plane  in  a  counterclockwise  sense  from  the  horizontal  direction  of 
strike.  These  angles  are  indicated  in  Figure  1. 

With  these  conventions,  all  possible  fault  planes  are  encompassed  by 
the  ranges  in  the  angles  of  0°  <  0t  <  360°,  0°  <  S{  <  90°,  and  -  180°  <  -If  <  180°. 
With  this  notation,  the  sense  of  P-wave  first  motion  at  the  center  of  the  focal 
sphere  is  positive  for  positive  values  of  Af  and  negative  for  negative  values. 
The  components  of  the  vectors  (Pujol  and  Herrmann,  1990)  are 
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Fig.  1.  Fault  plane  angle  convention.  The  x1(  xz  and  x3  axes  are  in  the  north,  east  and  down 
directions.  $  is  the  strike,  measured  north,  S  is  the  dip,  measured  downward  from  a  horizon¬ 
tal  direction  perpendicular  to  the  strike,  and  X  is  the  rake  angle  indicating  the  direction  of 
motion  on  the  fault,  given  by  the  vector  s.  The  side  of  the  fault  nearest  the  viewer  will  move 
in  the  a  direction. 


fl 

=  cos  Af  cos  <p{  + 

sin  Af 

cos  5( 

sin^f 

k 

=  cos  Ar  sin^f  - 

sin  Af 

cosS{ 

COS 

fa 

=  -  sin  Af  sin  Sr 

nl 

=  -  sin  <p{  sin  Sf 

n2 

-  cos0fsin£f 

n3 

=  -cos  S( 
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Another  way  to  combine  the  Green’s  functions  is  to  use  a  moment  tensor 
representation.  Jost  and  Herrmann  (1989)  related  the  Green’s  functions  in 
the  formulation  of  Her»-niann  and  Wang  (1985)  to  a  moment  tensor  source 
representation.  An  error  entered  into  the  Jost  and  Herrmann  (1989)  equa¬ 
tions  A5.4  -  A  5.6,  which  were  correct  only  for  a  pure  deviatoric  source  or  for  a 
pure  isotropic  source.  The  correct  expressions  are  as  follow: 


u£(r,  z,h,«y)  =  M** 

+  M 

+  M 

+  M 

+  M 

+  M 


ZSS  .  ZDD 
-  cos(2<j) - —  + 


y>: 


-ZSS  x  ZDD 
— =—  cos(2$>) - — 


ZDD  ZEP 
+ 


3  3 

ZSS  sin(2^)J 


xy 


y*i 


ZDS  cos(^)j 
ZDS  sin(^)j 


ZEP 

3 

ZEP 


(2a) 


ur(r,  z,h,6j) 


+  M 

+  M 

+  M 

+  M 


RSS  /o  ,  RDD  REP 
Tcot(2rt-  —  +-r 


cos(2«>)  - 


RDD 

6 


+ 


REP 

3 


RDD  REP 
3  +  3 


RSS  sin(2<>)  j 
RDS  cos(*)  j 


(2b) 
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+  My,  RDS  sin(^) 


u,(r,  z,  h,  o>)  = 


TSS 

2 


8in(2  0) 


+  Myy 


-TSS 

2 


sin(20) 


+  Mj-TSS  cos(2p)J  (2c) 

+  M„|tDS  sin(<>)j 
+  Mj-TDS  cos(*»]  . 

Here  the  moment  tensor  elements  are  with  respect  to  a  coordinate  system 
that  has  x-north,  y-east  and  z-down. 

The  displacements  corresponding  to  Green’s  functions  for  an  arbitrarily 
oriented  point  force,  given  by  the  vector  f  =  (flt  f2,f3),  are 

u,(r,  z,  h,  a)  =  (fj  cos  *  +  f2  sin  *)ZHF  +  f3ZVF 
Ujfr.z.h.aO^^cos^  +  f2  sin  0)RHF  +  f3RVF 

and 

u/r,  z,  h,  c o )  -  (ft  sin  £  -  f2  cos  #)THF  , 

where  the  1,  2  and  3  indices  refer  to  the  north,  east  and  down  directions,  and 
the  angle  </>  is  positive  from  north  to  and  east  direction.  The  meaning  of  the 
u„  ur  and  nt  is  the  same  as  above. 

Finally  the  PEP  solution  is  the  pressure  field  in  a  fluid  due  to  an  explo¬ 
sive  source  somewhere  in  the  model. 

WHOLESPACE  SOLUTION 

Explicit  expressions  for  the  Fj(k,  to)  functions  for  a  point  buried  source 
in  a  constant  velocity  wholespace  with  compressional  velocity,  a,  shear  veloc¬ 
ity,  0,  and  density,  p,  are  derived  from  Haskell  (1963,  1964)  as  follow: 

Define  the  vertical  wavenumbers  for  P-  and  S-  waves  as 
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{  Vk^kf  ksk„ 

Va  {iVkF^  k  <  ka 

and 

Vk2~k4  ksk, 

i^k|-k2  k  <  k,  ' 

For  this  case,  the  important  parameter  is  the  vertical  distance  between  the 
source  and  the  receiver.  Define  this  distance  tobeh  =  h-z.  A  negative  value 
of  h  indicates  that  the  receiver  is  beneath  the  source;  a  positive  value  indi¬ 
cates  that  the  source  is  beneath  the  receiver.  The  functions  appearing  in  the 
integrands  of  equation  (1)  are  as  follow: 


F1(k,e>)=  k  ,[(2k2  3k2)e-v»lhl  +  3k2e"v',hl]  sgn(h) 
AnpaP 

(3a) 

k  e-*'*  1  h  1 

F2(k,  o))=-  - - 5  [(2k2  -3k2) - +3v/>e-'' 1  h  *  ] 

A  xpa)1  va 

(3b) 

k2  1  h  1 

F3(k,ft>)  =  — — — -  [2vae~1'" ,h  1  -(2k2-k2) - ] 

Anpco1  p  Vp 

(3c) 

F4(k,a>)  =  ,  1  ,  [2k2e-‘'*,hI  -(2k2  k!)e-^,h']  sgn(h) 

4  npco2 

(3d) 

F6(k, a)  =  i^2  [e- 1  h  1  - 1  h  1  ]  sgn(h) 

(3e) 

F6(k, a)  =  j— — o  t~“  e”v“ 1  h  1  -Vpe~v>  ,h  1  ] 

4  Kp(Ol  va 

m 

F7(k,ffl)=4^2e’Ka'hl  8gn(h) 

(3g) 

F8(k,  to)  =  - — —  e~*'a  1  h  1 

Aitpalva 

(3h) 

F9(k,  <y)  =  — ^  [v„e-'- 1  h  1  -  —  e"^ 1  h  1 ) 

Ajtpco 1  vp 

(3i) 
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Fl0(k,<w)  =  - - F [e"v- ,h  1  -e-^,h,]8gn(h) 


Fn(k,<y): 


4  xpct)2 

-k2 
4  npeo2 


[e 


■  v-  I  h  1 


■e~v#lhl]  sgn(h) 


Fl2(k,a;)  =  TT7^  e','“ ' h '  -  ' h '  ] 


4  Jtp(D2  v0 


F13(k,  a)  = 
FM(k,e>)  = 
F16(k,  c o )  = 
F16(k,<y)  = 


e~v'lhl  sgn(h) 


4  itpffi 


An  pfivp 
-1 

4  np02vfi 
k 

4  na2v„ 


»-*> 1  h ' 


-Vfi  I  h  I 


Ihl 


(assumes  a  fluid  wholespace) 


(3j) 

(3k) 

(31) 

(3m) 

(3n) 

(3o) 

(3p) 


The  function  sgn(x)  is  defined  as 

-1  x<0 
8gn(x)  =  0  x=0 

1  x>0 

Following  Haskell  (1963),  it  is  possible  to  obtain  analytic  closed  form 
solutions  of  the  wholespace  Green’s  functions  corresponding  to  (la)  through 
(lj)  by  taking  partial  derivatives  of  the  Sommerfeld  integral  Fv: 


1  z™*  7 

F«  =  Re  ’  =1 


where 

and 


R 


R2  =  r2  +  h2 


e"‘,,lhlJ0(kr)dk 


v2  =  k2 


■6 


-9- 


(4a) 


The  closed  form  solutions  are  as  follow: 
1 


ZDD 


4  xpo)2  ah 


ah 


“  "Wi/,-3k,2dF^ 


ah" 


ah 


RDD  ~  A  1  2  [3  ^f0  3  +  ka2  9F°] 

4  xp<o2  ah2dr  ah23r  3r 

(4b) 

zds-  1  r°  ^F“  °  ^F/J  it  2^h 

4xp0)2  “  ah2dr  ”  ah2ar  *  5r 

(4c) 

f£DS  _  *  fo  ^  F“  o^F4»  ^  2  -I 

3rm  '  dh  ' 

(4d) 

X  2  a‘p,  a"F,  2aF, 

4 Kpo2  r  3r3h  ardh  p  3h  J 

_  i  a3Fa  a3Fa  2  aFa 

S  ■ttpai'1  ~  3r»3h  ah3  ”  ah 

(4») 

„  S3F,  c?Ffi  23F, 

ar2ah  ah3  '  ah1 

(40 

rss'4^[2  ar'  +ah!ar+k"  ft 

8%  ,3F, 

ar3  ah2ft  '  ft"1 

TSS  -  — (2  ^  +  2  ^  +  2k„2  ® “ 

4xp(o2  ar3  ah  3r  dr 

(4g) 

„S3Fa  „3*F,  2SF» 

ar3  “  ah2ar  ’  dr ] 

(4h) 

ZEP  -  -  1  aF“ 

4 npa2  ah 

(4i) 

REP  =  1  dFp 

4npa2  dr 

(4j) 

OT-  > 

4zpat2  dh  ah 

1  ,32F.  32F, 

L  "N  *M_  "V  -» 


RVF  = 


4 icpco2  drdh  drdh 
ZHF=-^[^f-^3 


4 npa>2  drdh  drdh 

RHF= _ ?h  k2Fi 

RHF  4np<o2 1  dr2  dr2  *  F*] 

1  1  aFa  dF,. 


THF  = 

PEP  = 


4/rp<u2  ^  r  ^  dr 

— - — F 
4  icpa>2  a 


dr 


•]  -k/F,] 


(4k) 

(41) 

(4m) 

(4n) 

(4o) 

(4p) 


where  the  partial  derivatives  are  evaluated  using  the  following  analytic 
expressions: 


dFv 

dh 

dFv 

dr 

d2Fv 

drdh 

d2Fv 

dh2 

d2Fv 

dr2 

d3Fv 

dr3 


-iiuR 

=  -  e  v 


-ia»R 


=  -  e  v 


-imR 


h  ie>  h 
R^+  v  R2 


r  iw  t 
R3+(V  R2 


=  -  e  v 


-  uuR 


f-Srh^  ,\co Y-Srh^j  ia>./-rh 


=  +  e  v 


-  icuR 


3h2  1 

R6  R3 


+  (- 


=  +  e  v 


-i«»R 


=  -  e  v 


(  9r  15rM  . 

[  R6+  R7  J+ 


/3h‘ 

1 

R4 

V 

~R2 

/ 3 r2 

1  ' 

w 

"  R2  J 

i<y/ 

9r 

_ 4, 

TV 

R4 

+  (-)21 
v 


V 


h2 

R3 


15r3 

R6 


-11- 


d3Fv 

ah3 


a3Fv 

3r23h 


a3Fv 

8h2dr 


3r  6r^ 
R3  +  R6 


-u»R 


e  v 


„  h  15h 
9R5+  R7 


3 


) 


ICO 

+  (  — 
V 


o  h  15h' 

9R4  +  R6 


3  ^ 


V 


V 


„  h  6h 
3  R3  +  R6 


•w-' 


J 


R4 

v  ) 


-3  Jl  +  15r*hV  ,i 


R5  '  R7  j"v\  ~  R4  +  R« 


-  iiuR 


e  v 


^  3r  15rh2)  ,i coi 
R®  +  ~ R7” 


v 


J 


3r  15rhz 
R4  +  R6 


(— )2I 

v 


'  r  6rh2  ^ 
R3+  R6 

v  J 


,ig\3 
+  (— )3 
v 


(  «  2  ^ 


rh4 

R4 


HALFSPACE  SOLUTION 

The  Fj(k,  co)  functions  for  a  point  buried  source  in  a  constant  velocity 
halfspace  with  compressional  velocity,  a,  shear  velocity,  p,  and  density,  p,  are 
derived  from  Haskell  (1963,  1964).  Let  h  and  z  be  the  depths  of  the  source 
and  receiver  beneath  the  free  surface,  respectively.  In  the  coordinate  system 
used,  both  are  positive  quantities.  Define  the  following  functions: 


Rayleigh  wave  period  equation: 
where 

y  =  2k2  /  k.0  . 

Free  surface  reflection  coefficients: 
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Rpp=[(r-D2+Aa^/k2]/FR 

RR  _  nZ 

pp  -  -Kpp 

R£s=-2r(r~l)/F* 

Rre  =[2r(r-lVa^/k2]/FR 

j»Z  _ pR 

"SP  -  lVj>g 
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F4(k,  eo) 

F5(k,o>) 

F6(k,<w) 

F7(k,/») 
F8(k  ,o)) 
F9(k,<w) 


-  — — y  (e-^lh-ll+I^se"^(h+1)+R|pe-('''h+,'*t) 
v„ 


)] 


4  Xp(D2 

2k2(e 1  h  -  * 1  sgn(h  -  z) + R^p  e  -  ^ ( h + 1 }  +  R^s  e  - ( *'-h + ^ } ) 
(2k2-k^Xe'‘'/,lh"z,8gn(h-z)+R^e”^(h+l)  +R|Pe''<v'h+‘''z))J 

4  npo)2 

(e-’'“lh-z,sgn(h-z)  +  RpP  e-v«(h+*)  +  RpSe~(v*h+*''z)) 
_(e-'.|h-*'sgn(h-z)  +  R|9e-,,'(h+z)+R|pe-(‘'«h+^z))] 


4xpco2 


[ 

k2 


.(g-.Jh-^l  +R«pe-Mh  +  z>+I^eMv0h+K,I)) 

u 

-  v/e~‘''lh“zl  +R|se_’''(h+z)  +  RgPe_(,,'h+‘',,z))] 

— (e-v“'h-zl8gn(h-z)+PZpe-,'»(h+z)+I^se-(,'»h+,''z)) 
4  rcpa2 

-k 

4  xpa2vt 


■  (e~ ,h -* 1  +  R*,  e-- (h+ z >  +  R?s  e ~ ( ’'•h+ v'l) ) 


\xp(o 2 


B(e“v'lh-Zl  +Rppe~v“(h+z)  +R|se'<v-h*^z)) 


k2 

_t_(e"'o,h-zl  +R|se-^<h+z)+R|pe-<v'h+l,"l) 


)] 


(5d) 


(5e) 


(50 


(5g) 

(5h) 

(5i) 
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(5j) 


F16(k,o)  =  - — - —  (e-''“lh~zl  -e_,'a(h+l))  (assumes  a  fluid  halfspace)  (5p) 
4nazva 

PLANE  LAYERED  MEDIA 

The  formulation  for  the  seismic  wave  field  for  a  buried  source  and  a 
buried  receiver  will  be  presented  for  a  medium  consisting  of  a  stack  of  layers 
within  two  halfspaces.  This  formulation  will  be  general  enough  to  encompass 
the  boundary  conditions  of  rigid  or  free  boundary  surfaces  and  also  elastic  or 
liquid  halfspaces. 

The  more  difficult  P-SV  problem  will  be  discussed  first,  with  a  simple 
extension  to  the  SH  problem  presented  next.  The  major  difference  between 


the  two  is  the  use  of  4x4  matrices  and  compound  matrices  rather  than  2x2 
matrices.  Within  each  section,  the  elastic  problem  will  be  discussed  first, 
together  with  a  simple  extension  to  a  liquid.  The  only  restriction  of  the  liquid 
model  with  respect  to  the  elastic  model,  is  that  a  liquid  layer  sandwiched 
between  two  solid  layers  will  not  be  permitted. 

Wang  and  Herrmann  (1980)  showed  that  the  Haskell  (1964)  formulation 
for  a  wavefield  due  to  a  buried  point  source  in  cylindrical  coordinates  could  be 
expressed  as 

®n-i  =  aN-i(dN_i)aN_2(dN_2) ' '  *  am(dm  ~  hm)AB  (6) 


+  aN-i(df<j_j  )aN-2^d>j_2)  •  •  •  ajld^Bo 

where  the  layering  convention  of  Haskell  ( 1964)  is  followed.  The  transformed 
motion-stress  vector  Bk  is  evaluated  at  the  k'th  interface.  For  P-SV  waves 
Bk  =  (Ur>  Uz,  Tz,  Tr)T  and  for  SH  waves  Bk  =  (U#,  T,)t  The  medium  properties  of 
the  k'th  layer  are  between  the  k  -  l'th  and  k'th  interfaces.  The  layer  thick¬ 
ness  of  the  k'th  layer  is  dk,  the  compressional  and  shear  wave  velocities  are 
«k  and  respectively,  and  the  density  is  pk.  The  source  is  at  a  depth  of  hra 
in  the  m'th  layer.  This  is  indicated  in  Figure  2.  Equation  (6)  states  that  the 
wavefield  at  some  position  beneath  the  source  is  a  function  of  the  wavefield 
discontinuity  at  the  source  and  also  the  wavefield  at  the  top  boundary.  There 
is  nothing  in  this  formulation  restricting  the  O'th  interface  to  be  a  free  sur¬ 
face,  or  the  N'th  medium  to  be  a  halfspace.  The  source  is  at  a  depth  hm  in  the 
m'th  layer. 

Haskell  (1964)  expressed  (1)  in  terms  of  upgoing  and  downgoing  P-  and 
SV-  potential  coefficients,  giving  a  total  of  four  coefficients  for  each  layer.  At 
the  top  and  bottom  boundaries,  the  boundary  conditions  require  only  two 
unknowns.  This  observation  is  used  below. 

P-SV  Problem 

Let  us  assume  that  there  are  matrices  G  and  H  such  that 


O' 

0 

a 


=  GBN1 


b 


and 


(7) 


B0  =  H(c,  d,  0, 0)T  (8) 

Multiplying  (6)  on  the  left  by  G  and  defining 
X  =  GaN_i(dN_i)aN_2(dN-2)  •  •  *  am(dm  -  hm) 
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Fig.  2.  Model  of  layered  medium,  showing  the  stress-displacement  vectors  at  the  interfaces 
and  the  medium  parameters  within  the  layers. 


and 

R  =  GaN_  i  (dfv_ !  )aN_2  ( dfj-2 )  • 
and  using  (7)  and  (8) 


■ai(di)H , 


O' 

c' 

0 

=  XS  +  R 

d 

a 

0 

,b_ 

0. 

(9) 


where  we  have  set  S  =  AB. 

Consider  the  first  two  equations  in  this  linear  system  of  four  equations. 
The  left  side  is  zero,  and  we  can  solve  directly  for  the  unknowns  c  and  d: 


(10) 


'C' 

-1 

R22  -R12  X..S, 

.d. 

^11^22  “1^12^21 

-R21  Rn  X^Si 
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~1  R22XliSi-R12X2iSi 

Rill  RaiXuSi-RuXaSi 


(11) 


~1  X2jZj2XliS1-XljZl2X2iSj 

Rill  "X^Z^X^Sj  +  XjjZjiX^Si 


(12) 


SiXl^Z^ 

-S|XlJ% 


(13) 


where  we  have  used  R  =  XZ,  and  the  compound  matrix  definition 
Rl{,  =  RikRji-RiiRjk.  Note  that  the  repeated  indices  represent  summations  in 
the  range  1-4. 

This  is  essentially  the  derivation  given  in  Wang  and  Herrmann  (1980), 
where  the  G  matrix  was  designed  to  give  the  upward  and  downward  poten¬ 
tials  in  the  halfspace,  and  the  H  matrix  was  such  to  make  the  free  surface 
stress  free,  with  c  =  Ur,  and  d  =  Uz,  the  displacements  at  the  top  boundary. 

The  formulation  of  (13)  has  proven  to  be  numerically  stable.  To  get  the 
wavefield  at  points  above  the  source,  it  is  obvious  through  the  use  of  propaga¬ 
tor  matrices  that  at  a  depth  hk  in  the  k'th  layer,  which  is  above  the  source. 


c 


““  )H] 


d 

0 

0 


(14) 


(15) 


One  may  be  tempted  to  use  (11)  first,  and  then  (12)  to  evaluate  the  wavefield 
of  a  buried  receiver  directly,  but  this  is  not  numerically  stable,  due  to  possibly 
increasing  exponential  terms  in  y  (In  the  extreme  case  of  a  layered 
wholespace,  the  whole  space  solutions  indicate  that  an  exponentially  decreas¬ 
ing  solution  is  required  in  the  z-direction  away  from  the  source). 

To  work  around  this  problem,  we  note  that  R  =  XZ  =  Xxy,  where  we 
define 


*  =  &m(hm)-  ■ak(dw-hk) 


(16) 


and 


y  =  ak(hk)"  a^d^H, 


(17) 
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where  me  propagator  matrix  property  a(h1)  a(h2)=a(hi  +  h2)  is  used.  Com¬ 
bining  equations  (13  -  17),  we  have,  after  some  algebra, 


Bfc  = 


IV 

Uz 

-1 

Tz 

11 

W 

10  fo 

TrJ 

L  L 

SiXl?,2  x,„  y! 
S,X|?,2  x,m  yl 
SiXl’2  X]m  yl 
S,Xli‘,2  xto  yl 


lm 

12 

2m 

12 

3m 

12 

4m 

12 


(18) 


If  this  is  written  in  terms  of  a  matrix  multiplication,  and  using  the  properties 
of  the  compound  matrix  that  Aljj,  =-A|Jk‘[  =-Al^  ,  and  that  Aljj,  =  0  for  i=j 
or  j  =  k, 


-1 


[Ur,  Uz,  Tz,  Tr]k  =  jfis  [Sl  S*.  S3,  S4] 


(19) 


0  Xl i|  X|  “  xl}2 
-Xlil  0  xlg  xlg 
-Xl  13  -X|23  0  Xl*2 

-XlJJ-Xl^-Xlg  0 


XU  *12  xia  x14‘ 

X21  *22  X 23  x24 
X31  X32  X33  X34 
-x41  X42  x43  x44. 


0  -ylM-ylu-ylu 
yl  il  0  -y I fi  -y 1 12 
yl  il  yl  ?!  0  -yIS 
yl  12  yl«  yl  12  o 


Noting  that  the  compound  matrix  indices  comprise  six  doublets,  we  can  fur¬ 
ther  simplify  the  expression  (14)  by  expressing  the  compound  matrix  ele¬ 
ments  associating  the  the  compound  matrix  doublets  {  12,  13,  14,  23,  24,  34  }, 
with  the  indices  {  1,  2,  3,  4,  5,  6  },  respectively,  defining,  for  example, 
X12  =X|  J|.  Thus  (14)  becomes 


[Ur,  U„  Tz,  Tr]k  =  ^|(Slt  S2,  S3,  S4) 


(20) 
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0 

Xu 

X12 

X13 

-Xu 

0 

Xu 

x15 

-X12 

-XM 

0 

X» 

~Xl3 

-Xj.5 

-X16 

0 

*11  *12  X13  X14 
X21  *22  X23  *24 
X31  X32  X33  X34 
X41  X42  X43  X44J 

•  0 
yn 
Y21 

.y3i 


•yn  -yai  -y3i 
0  -y4i  -ysi 
Y4i  0  -y61 
ysi  yei  0 


G  = 


Form  of  the  G  matrices  for  various  bottom  layer  conditions. 

Bottom  Halfspace  Free 

0  0  1  O' 

000  1 
1000 
.0  1  0  oj 

for  which  a  =  Ur,  b  =  U2,  and  we  force  Tz  =  0  and  Tr  =  0  at  the  N-l’st  interface. 
Bottom  Halfspace  Rigid 

'10  0  0' 


(21) 


G  = 


0  100 
00  10 
0  0  0  1 j 


(22) 


for  which  a  =  Tz,  b  =  Tr,  and  we  force  Uz  =  0  and  Ur  =  0  at  the  N-l’st  interface. 
Bottom  Halfspace  Elastic 


G  = 


~PY 

p(y-  D/va 

1  ~V}/va' 

1 

X- S, 

1 

b-* 

tT 

py!  k2 

Uvp  -1 

-py 

-p{y-  \)!va 

1  k2/Ka 

p{y-  D/vp 

py!  k2 

-If Vfi  -1  . 

(23) 


where  a  =  2A\  b=2B'  in  the  Haskell  (1964)  notation,  and  we  force  2A  "  =  0 
and  2B"  =  0.  This  guarantees  only  downwardly  propagating  wavefields  in  the 
halfspace. 
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Tb  form  the  compound  the  elements  of  the  R  and  X  we  need  only  the 
first  row  of  the  6x6  compound  G  matrix.  The  necessary  compound  matrix 
elements  of  the  G  matrix  are 


GIF 

Free 

Rigid 

Elastic 

GlJl 

0 

1 

p2^-r2/k2  +  (7-lf/ya^ 

Glif 

0 

0 

-plvp 

Glu 

0 

0 

p(r  -  kfy-D/y^j 

Gig 

0 

0 

p^-y  IV?  +  {y  -  \) /  vav p 

Gl24 

0 

0 

P>va 

GIL2 

1 

0 

k2/ vavp ~ 1 

Form  of  the  H  matrices  at  the  top  interface 

Top  Surface  Free 


H 


'10  0  0- 
0  100 
00  10 
0  0  0  1. 


(24) 


for  which  c  =  Ur,  d  =  Uz,  and  we  force  Tz  =  0  and  Tr  =  0  at  the  0’th  interface. 


7b p  Surface  Rigid 


H  = 


0  0  1  O' 
000  1 
1000 
0  100 


(25) 


for  which  c  =  Tz,  d  =  Tr,  and  we  force  U2  =  0  and  Ur  =  0  at  the  0’th  interface. 
Elastic  Medium  above  top  surface 
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2H  = 


(26) 


-1//7  Vp/(>  -1/p  -Vp/p' 

-va! p  k %/p  va! p  k2/ p 
-<y~  l)  yvp  — (y — D  -yvp 
_-yva/k2  y- 1  yvj k2  y- 1  . 

for  which  c  =  2A",  d  =  2B",  and  we  force  2A'  =  0  and  2B'  =  0  at  the  O’th  inter¬ 
face.  This  guarantees  only  upwardly  propagating  wavefields. 

The  necessary  compound  matrix  elements  of  the  H  matrix  are 


Hii2 

Free 

Rigid 

Elastic 

Hi  if 

1 

0 

(vavfi-k2)/4p2 

HIS 

0 

0 

-Vp/4  p 

Hi  I, 

0 

0  ( 

yvaVp/k} -{y  -  l)j/4p 

Hl§ 

0 

0 

\)s.Hy-l)-yvav^/4p 

Hlg 

0 

0 

vaJ4p 

HlfJ 

0 

1 

y2vavp/k2 -{y  -lf^j4 

LIQUID  LAYER 

The  inclusion  of  a  liquid  layer  in  the  model  complicates  the  formulation, 
since  only  P-potentials  are  required  and  the  Haskell  matrix  formulation  in 
terms  of  P-  and  SV-arrivals  does  not  simply  reduce  to  the  fluid  problem.  The 
4x4  propagator  matrices  in  the  elastic  medium  must  then  be  connect  to  2x2 
propagator  matrices  in  the  fluid.  This  complicates  the  structure  of  the  pro¬ 
gram. 

If  a  simplifying  assumption  is  made,  that  the  liquid  layer  is  not  between 
two  elastic  layers,  the  mathematics  is  simplified.  This  restriction,  means  that 
the  4x4  matrix  computations  will  be  maintained,  but  with  a  slightly  different 
meaning  for  the  propagator  matrix.  Basically,  in  a  fluid,  only  U,  and  T,  must 
be  computed,  since  only  these  two  quantities  are  continuous  at  boundaries.  In 
a  fluid  Tr  is  by  definition  zero,  and  Ur  is  discontinuous  at  boundaries. 
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Let  the  1  +  l'st  layer  be  elastic  and  the  l'th  layer  be  fluid.  Let  the  elastic 
displacements  at  the  top  of  the  elastic  layer  be  related  to  the  fluid  displace¬ 
ments  at  the  top  of  the  fluid  layer  by  the  following  pseudo-propagator  matrix 
relation: 

-uri  r  i  o  o  o]r  air\  - 

U2  _  0  cosh  vad  ~va  sinhvad/p  0  (U,),.! 

Tz  0 -psinhvad/va  coshvad  0  0^)^ 

TrJ,  Lo  0  0  lJL  (Tr),  _ 

This  pseudo-propagator  propagates  the  vertical  displacment  and  stress 
downward  to  the  top  of  the  elastic  layer  for  future  use,  while  placing  no  con¬ 
straint  upon  the  radial  displacement  or  stress.  This  was  first  used  by  Hud¬ 
son  (1969).  For  example,  the  requirement  that  (Tr)j  =  0  can  be  met  after  mul¬ 
tiplying  by  the  4x4  matrix  in  (27).  The  restriction  that  the  fluid  layer  cannot 
be  between  two  solid  layers,  arises  from  the  fact  that  four  free  parameters, 
radial  displacement  and  stress  at  the  two  elastic  fluid  boundaries,  must  be 
saved,  but  this  formalism  permits  only  two  free  parameters  to  be  retained. 
Tb  compute  Ur  in  the  fluid,  we  note  that  it  is  proportional  to  Tz  in  the  fluid. 

Given  this  preface,  the  following  matrices  are  used  to  meet  the  possible 
boundary  conditions  if  a  fluid  layer  is  at  the  surface  or  the  lower  halfspace  of 
the  model: 

The  G  matrix  for  the  halfspace. 

Bottom  Halfspace  Free 

00  10 
0001 
1000 
0  100 

for  which  a  =  Ur,  b  =  Uz,  and  we  force  Tz  =  0  and  Tr  =  0  at  the  N-l’st  interface. 
Bottom  Halfspace  Fluid 

0  0  0  1 

0  -p/2va  1/2  0 
0  p/2v„  1/2  0 
10  0  0 

where  a  =  A',  b  =  Ur,  and  we  force  A"  =  0  and  Tr  =  0.  This  guarantees  only 
downwardly  propagating  wavefields  in  the  halfspace. 

The  necessary  compound  matrix  elements  of  the  G  matrix  are 
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Gig 

0 

0 

Gill 
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0 

GIL2 

0 

0 

Gl  23 
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0 

GIL2 

0 

p!  2va 

GIL2 

1 
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Form  of  the  H  matrix  at  the  top  interface 
Tbp  Surface  Free 


H  = 


1  0  0  O' 
0  100 
0010 
0  0  0  1. 


(30) 


for  which  c  =  Ur  at  top  of  closest  elastic  layer  beneath  the  fluid,  d  =  U,  at  the 
top  of  the  fluid,  and  we  force  T,  =  0  at  the  top  of  the  fluid  and  Tr  =  0  at  the 
top  of  the  closest  elastic  layer  underneath.  The  top  of  the  fluid  is  the  0’th 
interface. 

Fluid  halfspace  above  top  layer 


r  1  0  0  01 


2H  = 


0  -ya / p  va/p  0 
0  1  10 

0  0  0  1. 


(31) 


for  which  c  =  Ur  at  the  closest  elastic  interface,  d  =  A",  and  we  force  A'  =  0  at 
the  0’th  interface  and  Tr  =  0  at  the  top  of  the  nearest  elastic  boundary.  This 
guarantees  only  upwardly  propagating  wavefields  into  the  upper  halfspace. 

The  necessary  compound  matrix  elements  of  the  H  matrix  are 


Hi  i2  Free  Elastic 
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HIS  1  -vjp 

Hlg  0  1 

Hi  g  0  0 

Hlg  0  0 

Hlg  0  0 

Hlg  0  0 


SH  Problem 

The  SH  development  uses  the  same  formalism  as  (6),  except  that  the  ak 
matrices  are  now  2x2  rather  than  4x4.  To  represent  the  boundary  conditions 
at  the  top  and  bottom  boundaries,  we  assume  that  there  are  g  and  h  matrices 
such  that 

°  =gBN  (32) 

e 

and 


(33) 


Here  B  =  ,T,J  .  Multiplying  (32)  on  the  left  by  G  and  using  (6),  and  defin¬ 
ing 

x  -  fSas- i(dN-i )aN_2(d}sr_2)  •  •  •  am(dm  - hm) 


and 


we  have 


r  =  gaN_i(dN_1)aN_2(dN_2)  •  •  ■  a^d^h  , 


where  we  have  set  s  =  AB. 


0 

e 


=  xs  +  II 


f 

0 


(34) 
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The  coefficient  f  is  quickly  found  by  the  relation 

f_  |*ngi +  *i2S2 j 


(35) 


Tb  find  the  displacement  field  anywhere  between  the  top  boundary  and  the 
source  depth,  we  use  the  propagator  matrix 


Bk  =  afc(hk)---ai(di)] 


=y 


Thus  the  transformed  tangential  displacement  is  just 

/XnS1+Xa2S2'| 

u*='y“l  TTi  J’ 

which  is  a  stable  numerical  function. 

Form  of  the  g  matrix  for  various  bottom  layer  conditions 

Bottom  Boundary  Free  Surface 

0  i-1 


(36) 


(37) 


(38) 


*  = 


1  0 


(39) 


Bottom  Elastic  Halfspace 


g  = 


PVfl 

pVfi 


_1_ 

-1 

p  J 


(40) 


(The  true  expression  should  have  terms  like  to2 1  ft 2  instead  of  1/fi2,  but  the 

source  expressions  are  either  of  the  form  a  0]T  orjo  b/<y2j  .  Thus  the  multi¬ 
plications  required  for  (38)  will  eliminate  the  0?  terms,  and  the  apparent 
o>  =  0  singularity  is  numerically  avoided  by  using  this  form  for  the  g  matrix 

and  [0  bj  for  the  source  expression.) 

Bottom  Halfspace  Rigid 


1  0 
0  1 


(41) 
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Form  of  the  h  matrix  for  various  top  layer  conditions 

Top  surface  free 


Tbp  surface  elastic 


Tbp  surface  rigid 


0 

1 


1 


1 

0 


(42) 


(43) 


(44) 


Source  Terms 

The  source  terms,  S,  representing  the  discontinuity  in  the  displace¬ 
ment-stress  values  across  the  source  layer  are  given  in  the  following  table  for 
different  source  representations.  The  S*  values  are  the  P-SV  source  terms 
used  in  equations  (10-13,  18,  19)  while  the  Sj  are  the  SH  source  terms  given 
in  Table  1. 
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Term 

Si 

Table  1. 

S* 

Source  Term  Coefficients 

Sa  S4 

8l 

82 

DD 

0 

2ki 

4kk2 
Ax  pat2 

0 

2k^(20/a)2-3j 
Ax  to2 

0 

-2 

0 

DS 

P 

Axpofi 

0 

0 

0 

-2k 

Axp/3 2 

0 

2k 

SS 

0 

0 

2kk2 

0 

Axo)2 

4kk2 

0 

4* 

EP 

0 

Ax  pa2 

0 

-2k 

4*<»2kJ 

0 

0 

VF 

0 

0 

Axo)2 

0 

-2 

0 

0 

2 

HF 

0 

0 

0 

Axo)2 

0 

Ax 

lb  obtain  the  required  integrand  in  (1),  the  above  source  coefficients  are 
used  with  (20)  for  Uz  and  Ur  and  (38)  for  U#  to  yield  the  following  Fj  factors. 
By  taking  the  negative  of  the  U*  terms,  we  guarantee  that  the  corresponding 
Green’s  functions  yield  positive  ground  motion  in  an  upward  direction. 
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Source  Term 


N 
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Sdd 
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Sdd 

f3  =  -uz 

Sds 

f4  =  ur 

Sds 

f«— u. 

Sss 

f6=  ur 

Sss 

F7  =-Ut 

Sgp 

Fa=  Ur 

Sep 

Fg  =  -Ut 

Svf 

F  io  =  ur 

Svf 

Fn=-Uz 

Shf 

f12  =  ur 

Shf 

f13  =  u, 

SDS 

F*  =  U, 

8ss 

Fis  =  U, 

SHF 

Fie  =  pF8 

(receiver  in  fluid) 

Evaluation  for  receiver  beneath  the  source. 

The  matrices  used  for  the  general  solution,  (20)  and  (38)  were  derived 
for  the  source  position  beneath  the  receiver.  It  was  because  of  this  assump¬ 
tion  that  the  relatively  simple  expression  were  obtained.  If  the  source  is 
above  the  receiver,  then  an  alternative  derivation  is  required,  first  to  obtain 
the  coefficients  c,d  or  e.  From  these,  Bn_x  is  found,  and  then  the  propagator 
matrices  can  be  used.  The  problem  with  this  technique  is  that  extra  compu¬ 
tational  steps  must  be  placed  into  the  coding.  The  alternative  approach  is  to 
invert  the  order  of  the  layered  model  by  reversing  the  sense  of  the  z  axis,  to 
note  the  source  and  receiver  positions  in  the  new  model,  and  then  to  evaluate 
the  solution  using  (20)  and  (38).  Finally,  care  is  taken  to  preserve  the  sense 
of  positive  positive  displacements  of  the  original  model.  This  Iasi  ask  is 
accomplished  in  two  stages.  First  the  S2,  S4  and  s2  coefficients  are  rep*aced  by 
~S2,  -S4,  and  -s2.  since  these  are  the  source  coefficients  that  depend  on  the 
direction  of  the  local  vertical.  Once  the  final  displacments  are  computed,  the 
vertical  displacements  are  inverted,  while  the  radial  and  transverse  time  his¬ 
tories  are  unchanged.  This  sequence  of  steps  is  tested  by  looking  the  the  P- 
and  S-wave  particle  motions  for  receivers  above  and  beneath  the  source. 

TIME  HISTORIES 

Evaluation  of  the  inverse  Fourier  tranfrom  to  yield  a  time  series  is  usu¬ 
ally  accomplished  by  an  inverse  Fast  Fourier  Transform  (Brigham,  1974) 
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which  approximates  the  true  inverse  transform 

oo 

g(t)=  J  G(f)ei2*ftdf 

—  oo 

by  the  inverse  Discrete  Fourier  Transform 

g(kAt)  =  £ GinADe^^^Af  for  k  =  0, ....  N- 1 

n  =  0 

where  Af=  (Note  that  this  is  related  to  the  original  Fourier  transform 
NAt 

definition  defined  in  the  introduction  if  h(t)=g(t)  and  G(f)=H(2*f)).  This 
approximation  to  the  continuous  Fourier  Transform  introduces  concerns 
about  finite  frequency  and  time  windows  and  periodicity  in  the  time  and  fre- 
quenccy  domains.  This  will  always  be  present  when  the  discrete  Fourier 
transform  pair  is  used,  but  can  be  controlled  so  that  the  discrete  Fourier 
transform  result  is  a  reasonable  approximation  to  the  desired  Fourier  trans¬ 
form  solution. 

The  first  problem  considered  concerns  the  effect  of  the  finite  frequency 
window  in  the  resultant  time  histories.  The  function  H(f)  is  sampled  in  the 

range  (-fN,fN)  where  =  Frequencies  outside  this  range  are  effectively 

set  to  zero.  If  the  transitions  through  the  frequency  points  f  =  ±fN  are  not 
smooth,  high  frequency  ripples  will  be  seen  in  the  time  domain.  This  effect 
can  be  significantly  reduced  in  synthetic  seismograms  by  choosing  source 
time  functions  that  have  zeros  at  the  Nyquist  frequencies  ±fN.  Two  possible 
functions  that  have  this  property  are  the  triangular  and  parabolic  pulses. 


Triangular: 


s(t)  =  -  ( 
X 


which  has  the  Fourier  transform 


0 

t 

T 

1-1 


0 


t  <0 
0<t  <r 
r  <  t  £  2r 
t  >  2r 


sin  *fr 
/rfr 

L 

This  is  a  positive  pulse  with  unit  area  and  a  comer  frequency  4  =  l/*r.  If 
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r  =  2MAt  and  M>  1,  this  function  has  spectra]  zeros  at  frequencies  ~  fN, 

2  1  M 
—  ffi,  ....  fN,  where  fjj  is  the  Nyquist  frequency  defined  as  fN  =  Thus  if 

we  sample  the  Fourier  transform  of  this  pulse  and  apply  the  inverse  discrete 
Fourier  Transform,  we  will  see  a  nice  sampled  triangular  pulse  without  any 
ripples  in  the  time  domain. 


Parabolic: 


s(t)  = 


2r 


0  t<0 

l/2(t/r)2  0<t<r 

-  l/2(t/r)2  +  2(t/r)-  1  r<t<3r 
l/2(t/r)2-4(t/r)  +  8  3r<t<4r 

0  t  >  4r 


.  2  0)7 

4  sin  —  sin  cdt 


(cot)3 


The  Fourier  transform  of  this  function  is 

r-ia>2t  [2  sinair-sin  2o)t]  __  ia)2r 
(air)3 

This  time  function  has  a  unit  area  and  a  comer  frequency  fc  =  1/4.575 r.  In 
addition,  it  has  spectral  zeros  at  certain  frequencies.  If  r  =  MAt,  where  M  >  1, 

1  2 

then  spectral  zeros  are  at  frequencies  —  fN,  —  fN,  . . . ,  fN.  By  choosing  r  and 

M  M 

At  such  that  one  of  the  spectral  zeros  occurs  at  the  Nyquist  frequency,  the 
pulses  can  be  synthesized  and  propagated  through  the  model  without  the  rip¬ 
pling  introduced  by  an  arbitrary,  sharp  high  frequency  spectral  cutoff.  Note 
that  the  parabolic  pulse  with  r- At  will  give  the  same  sample  values  at  the 
triangular  pulse  with  r=2At. 

The  other  problem  to  be  addressed  is  that  of  the  periodicity  of  the  dis¬ 
crete  Fourier  transform  pair.  This  can  be  addressed  by  introducing  complex 
frequencies. 

Consider  the  the  Fourier  transform  pair 

oo  oo 

G (co)  =  [  gftJe-^dt  g(t)  =  ~  f  G(r,-)ei‘atd(« 


or  symbolically 


g(t)  <=>  G(<y) 


From  this  definition,  we  see  that 

e~aigit) 


G{(o-ia) 
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This  is  equivalent  to  stating  that 

1  00 

git)  =  eal  —  [  G(<y-ia)e’*'ld<2> 
2jt  J 

-00 

The  convolution  theorem  states  that  given 

x(t)  <=a  X(<a) 


then 


and 


y(t)  <=*>  Y (a>) 

z(t)  =  x(t)  *  y(t)  <=>X(a>)Y(at) 


e  z(t)  = 


e  x(t)  I  *  I  e  y(t) 


X(<u-ia)Y(<y-ier) 


The  time  shift  theorem  states  that 

g(t-to)  «=>  G(<y)e"1,ut<) 

and  hence 


e-a(t-io)g(t_t0)  <*=*>  G^-iale-’"10 
or,  by  using  the  linearlity  of  the  Fourier  transform, 

e'atgU-to)  <=*>  G((o-ia)e'(“-'a)u> 

This  indicates  how  a  time  shift  is  formed  in  the  frequency  domain  and  the 
proper  de-attenuation  factor  in  the  time  time. 

The  value  of  using  the  time  domain  attenuation  is  obvious  when 
attempting  to  propagate  a  simple  pulse,  e.g.,  to  evaluate  the  inverse  trans- 

.  <a 

form  of  H(<»)e  'cx  as  x  increases.  As  x  increases,  the  pulse  will  arrive  at  later 
times.  However,  it  will  eventually  wrap  around  to  appear  at  zero  time.  This 
periodicity  effect  can  be  reduced  by  using  the  time  domain  attenuation  tech¬ 
nique.  Upon  wrap  around,  the  pulse  will  be  attenuated  in  amplitude  by 
e-«NAt  obviously  the  larger  the  value  of  a,  the  lower  the  amplitude  of  this 
effect. 

An  important  added  advantage  of  this  technique  is  that  it  will  remove 
the  surface  wave  poles  and  branch  point  singularities  from  the  real 
wavenumber  axis  for  a  perfectly  elastic  problem,  so  that  contour  integration 
is  not  required  in  the  wavenumber  domain  (Phinney,  1965). 

An  example  of  this  technique  is  given  in  Figure  3.  For  computational 
speed,  the  RDD  Green’s  function  for  a  wholespace  (4b)  is  evaluated.  The  P- 
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wave  and  S-wave  velocities  are  6.0  and  3.5  km/s,  respectively,  the  density  is 
2.8  gm/cm3,  the  receiver  is  a  height  5  km  above  the  source,  and  synthetics  are 
computed  to  distances  of  300  km.  All  traces  start  at  -5.0  and  terminate  at  58 
seconds.  The  sampling  interval  is  1.0  s,  and  the  triangular  pulse  has  a  base  of 
4.0  s,  e.g.,  r  =  2. 0s.  Figure  3a  shows  the  result  of  using  a  -  0.0.  In  this  case 
the  wrap  around  of  the  S-wave  arrival  beyond  175  km  is  obvious.  By  using  a 
=  0.04  in  Figure  3b,  the  amplitude  of  the  wrapped  S-wave  is  reduced.  Of 
course,  if  the  P-wave  wraps  around,  the  traces  will  look  different,  and  the  P- 
wave  amplitude  will  be  significantly  reduced. 


Fig.  3.  Demonstration  of  the  effect  of  the  time  domain  damping  factor,  or  (a)  a  =  0.0  ;  (bj  a  = 
0.04.  The  RDD  Green’s  function  for  a  wholespace  with  Vp  =  6.0  km/s,  Vg  =  3.5  km/s  and  p  = 
2.8  gm/cm  .  The  velocity  traces  (peak  velocity  in  cm/s  is  shown)  show  the  Green’s  functions 
at  distances  of  0  to  300  km  for  a  source  buried  at  5  km  beneath  the  receiver.  All  traces  start 
at  -5.0  seconds  and  end  at  58.0  seconds.  A  triangular  pulse  with  base  of  4.0  seconds  is  used 
as  the  source  vavelet. 


The  choice  of  a  is  not  simple.  First,  if  a  is  not  large  enough,  then  the 
wrapped  signals  will  not  be  significantly  reduced  in  amplitude.  On  the  other 
hand,  if  a  is  too  large,  then  the  wavenumber  integrands  in  Figure  1  may  be 
so  smoothed,  that  numerical  integration  noise  may  become  important.  For 
models  without  significant  reverberations,  e.g.,  thin  surface  water  layers, 
reducing  the  wrapped  amplitude  by  a  factor  of  10  will  suffice.  Thus  if  the  time 
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window  is  NAt  seconds  long,  choose  a  such  that  aNAt  =  2.3.  In  the  example 
just  run,  this  product  was  2.52. 

WAVENUMBER  INTEGRATION 

The  evaluation  of  the  wavenumber  integrals  in  (1)  is  complicated 
numerically.  Numerical  evaluation  introduces  error  through  the  finite  inte¬ 
gration  limits  and  through  the  particular  integration  rule  used.  These  errors 
may  appear  in  the  time-distance  domain  as  propagating  arrivals  or  as  noise 
unrelated  to  distance.  Items  of  concern  are  the  infinite  limit  of  integration  in 
wavenumber,  the  integration  scheme  used,  and  phase  velocity  filtering. 

Infinite  Limits 

The  problem  is  to  approximate  the  infinite  integral  by  a  finite  one.  The 
obvious  approach  is  to  truncate  the  integral  so  that 

00 

Jflk,r)dk=  J  ftk,  r)dk, 

c  o 


but  now  the  choice  of  kmax  becomes  important,  especially  when  the  function 
flk,r)  is  significantly  different  from  zero  for  k  >  kmax. 

To  illustrate  this,  consider  the  REP  function  for  an  elastic  wholespace 
(lj,  3h,  and  4j) 


4npa2  dr 
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v«lhl  Jj(kr)kdk 


The  consideration  in  choosing  kmax  is  it  must  be  large  enough  to  give  the 
desired  signals.  If  the  I  h  I  is  not  zero,  then  there  will  be  some  value  of  kmax 
that  together  with  the  I  hi  will  reduce  the  integrand  significantly  so  that  the 
rest  of  the  integral  from  k  =  kmax  to  k  =  oo  can  be  sucessfully  ignored. 

Unfortunately,  this  cannot  be  done  practically  when  I  hi  is  small.  To 
handle  this,  one  can  make  use  of  the  asymptotic  value  of  the  integrand. 
Returning  to  the  expression  for  REP,  note  that  the  for  large  values  of  the 
wavenumber,  that  va  -  k  and  that  the  integral 
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This  is  just  the  expression  for  REP  evaluated  at  co  =  0.  Note  that  for  large  val¬ 
ues  of  the  wavenumber,  both  integrands  behave  similarly.  This  suggests  the 

oo 

following.  If  g(k,  r)  =  f(k,  r)  as  k  oo,  and  if  J  g(k,  r)dk  =  CXr),  then 

o 

oo  oo 

J  fik,  r)dk  =  j  £f(k,  r)  -  g(k,  r) jdk  +  G(r) 
o  o 


J  flk,  r)dk  =  J 

o  o 

where  k^  is  such  that  now  the  integral 


flk,  r)-g(k,  r)  dk  +  G(r) 


is  negligible. 

There  are  now  two  choices  to  be  made:  the  value  of  k^  and  whether  or 
not  to  use  the  asymptotic  integration  trick. 

The  choice  of  kmax  requires  a  tradeoff  between  being  accurate  and  being 
efficient.  At  low  frequencies,  the  choice  should  depend  on  the  vertical  distance 
term,  h.  On  the  other  hand,  at  high  frequencies,  the  real  part  of  vjhl 
becomes  large,  and  a  frequency  -  depth  dependent  limit  may  be  appropriate. 

One  strategy  that  seems  to  work  and  to  be  efficient  is  as  follows:  Define 
I  hi  to  be  the  absolute  value  of  the  difference  m  the  source  and  receiver 
depths  in  the  model.  Also  define  H  to  be  a  mean  thickness  of  layers  in  the 
model.  Also  let  kv  be  the  wavenumber  associated  with  the  minimum  wave 

*mm 

velocity  (usually  the  S-wave  velocity,  but  may  be  the  P-wave  velocity  of  a  fluid 
layer)  at  the  current  angular  frequency  cd  and  let  k^m  n  be  the  wavenumber 
associated  with  the  minimum  wave  velocity  but  at  the  maximum  desired  fre¬ 
quency  in  the  synthetics,  co^. 


a) 

b) 


The  steps  are  now  as  follow. 

5 

Estimate  a  test  variable  k=  —  +  FAC  kv  . 

m  U  vmin 


if  kmh  >  5 


,  6.0 

1c  — 
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else 


kj  =  20k? 
k2  =  5k? 
ifkjhl  >5 
,  6.0 
ki  = 


k2 

km, 


Ihl 
=  k2 


+  4k„ 


if  kv  >  kj 

vmtn  x 

Do  not  use  asymptotic  integration  technique 

else 

Do  use  asymptotic  integration  technique 

The  purpose  of  this  logic  is  simple.  If  the  user  specified,  frequency 
dependent  upper  limit,  km,  is  large  enough  then  the  expC-Vy  I  h  I )  term  is  small 
enough  so  that  the  truncation  error  is  negligible.  This  technique  is  described 
in  Apsel  and  Luco  (1983). 

If  this  is  not  true,  then  use  the  asymptotic  values  to  truncate  the  inte¬ 
gral.  Depending  upon  the  particular  Green’s  function,  the  function  g(k,  r)  is 
defined  to  be  of  the  form  (A+Bk)exp(-k|h|)  or  (Bk+Ck2)exp(-k|hl ).  The  two 
constants  are  estimated  using  the  ki  and  k2  values.  The  numbers  6.0  and  2.5 
were  chosen  since  exp(-6.0)  and  exp(-2.5)  differ  sufficiently  to  define  the  con¬ 
stants  but  neither  is  so  small  that  one  runs  into  lack  of  significance  in  using 
single  precision  arithmetic.  The  particular  form  chosen  for  km  is  designed  to 
be  economical  at  high  frequency,  but  also  to  guarantee  sufficient  sampling  at 
low  frequencies. 

When  the  difference  in  the  source  and  receiver  depths  is  small,  using  an 
upper  limit  proportional  to  Ih!"1  would  lead  to  excessive  computational  effort. 
Thus  a  decision  is  made  to  tie  the  computations  to  the  wavenumber  of  the 
highest  frequency. 

The  FAC  parameter  also  controls  the  upper  limit,  and  requires  some 

(  Y/2 

judgment.  If  FAC  is  made  too  small,  then  k„  -  k2  will  not  be  close  enough 

to  km,  in  which  case  a  truncation  error  will  introduce  spurious  low  velocity 
arrivals.  On  the  other  hand,  FAC>  1  is  required  to  include  a  Rayleigh-wave 
pole,  a  FAC  =  3  or  4  seems  to  work  well. 

If  one  considers  the  range  of  integration  in  the  two-dimensional 
omega -k  domain,  a  triangular  region  is  sampled  when  I  hi  is  large,  and  a 
somewhat  rectangular  region  when  it  is  small.  Since  the  upper  wavenumebr 
limit  depends  on  the  highest  frequency,  increasing  the  frequency  content  of  a 
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signal,  while  keeping  the  length  of  the  time  series  and  the  wavenumber  sam¬ 
pling  fixed,  will  increase  the  computational  effort  quadratically. 

Bouchon  Integration  Scheme 

Bouchon  (1981)  analyzed  the  use  of  a  trapezoidal  integration  rule  to 
evaluate  the  integrals  in  (1).  Bouchon  proved  that 


+ 


dk 


where  Ak=2;r/L,  r0  =  -  ,  £j  =  1  for  j  >  0,  and  kn  =  nAk.  By  considering  the  Som- 

merfeld  integral,  it  can  be  shown  by  large  wavenumber  approximations  to  the 
integrals  that  the  cosine  term  introduces  a  sequence  of  noise  arrivals  corre¬ 
sponding  to  rings  of  sources  at  distances  L,  2L,  which  introduce  waves 
propagating  toward  and  away  from  the  origin.  From  this  consideration,  Bou¬ 
chon  (1981)  recommended  that  the  parameter,  L,  be  chosen  such  that 


L>2r 


and 

r  i1/2 

l(L-r)2  +  z2J  >vt, 

where  z  is  the  vertical  distance  between  the  source  and  the  receiver,  r  is  the 
radial  distance,  t  is  the  maximum  time  for  which  the  trace  is  to  be  generated, 
and  v  is  the  velocity  of  the  wave.  When  these  rules  are  evaluated  for  a  more 
complicated  model,  then  the  z  represents  the  total  vertical  path  of  the  last 
significant  arrival  and  v  is  the  fastest  velocity.  The  choice  of  L  is  a  matter  of 
experiment. 

Herrmann  and  Mandal  (1986)  modified  the  Bouchon  (1981)  expression 
by  using  a  shifted  rectangular  rule  because  of  the  appearance  of  non-causal 
arrivals  (  a  k  =  0  contribution)  in  some  of  the  integrands, such  as  RDS.  This 
noise  arrival  may  be  inherent  to  the  use  of  a  rectangular  or  trapezoidal  rule, 
but  may  also  be  linked  to  the  scheme  used  to  truncating  the  limits  of  integra¬ 
tion.  They  used 
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+ 


2  cos  (n(k  -  ko)L))  )  dk 


where  Ak  =  2*7L,  r0=  fj  =  1  for  j  >0,  and  kn  =  nAk  +  k<).  The  value  for  ko  was 

z 

determined  empirically  to  be  ko  =  0. 218Ak. 

Other  Integration  Schemes 

The  use  of  the  Bouchon  (1981)  integration  scheme  is  not  computation¬ 
ally  efficient  at  large  distances  and  high  frequencies,  because  of  the  required 
sampling  in  the  k -a)  domain.  Another  approach  is  to  approximate  the  inte¬ 
grand  by  a  higher  order  polynomial  in  the  interval  from  k  to  k  +  Ak.  If  one 
part  of  the  integrand  is  oscillatory,  e.g.,  the  Bessel  function,  then  a  Filon  inte¬ 
gration  scheme  can  be  developed.  Apsel  and  Luco  (1983)  used  a  fourth  order 
polynomial  to  model  the  non-oscillatory  part  of  the  integrand. 

Mallick  and  Frazer  (1987)  suggested  the  use  of  Filon  rule  based  on  a  lin¬ 
ear  approximation  to  the  non-oscillatory  part  of  the  integrand.  Thus  if  the 
integrand  is  of  the  form 

k2 

J  F(k)e"ikrdk , 

k, 

it  can  be  approximated  in  the  case  of  r  5*  0  as 

where  <JX=X(k2)-X(k1).  Mallick  and  Frazer  (1987)  and  Saikia  (1993)  used 
this  form  together  with  the  asymptotic  expansion  of  the  Bessel  or  Hankel 
functions  to  evaluate  the  wavenumber  integrals  of  the  type  discussed  in  this 
paper. 

Integrands  for  Large  Offset 

Bouchon’s  (1981)  mathematical  development  demonstrated  the  there 
can  be  inwardly  propagating  noise  arrivals  because  of  integrating  the  Bessel 
function.  At  large  distances,  r,  cleaner  seismograms  may  arise  by  using  the 
following  approximation: 

oo  1  00 

J  F(k)Jm(kr)dk  =  -  J  F(k)H^(kr)dk, 
o  o+ 

where  (z)  is  the  Hankel  function  of  the  second  kind.  This  is  a  standard 
approximation,  which  can  be  justified  from  physical  grounds  that  only  out¬ 
wardly  propagating  signals  are  desired.  The  mathematical  justification  arises 
from  applying  both  the  physical  requirement  and  the  stationary  phase 
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approximation  to  drop  the  first  order  Hankel  function  from  the  identity 

2JJz)=H<»(z)+H<2)(z). 

Phase  Velocity  Filtering 

One  final  modification  of  the  integrand  is  due  to  Fuchs  and  Muller 
(1971).  This  is  the  concept  of  computing  synthetics  with  a  range  of  phase 
velocities  bounded  by  [Cm^Ci.c^CmiJ.  These  can  be  used  to  define  a  window 
function,  W(k)  defined  by 


W(k)  = 


1-cos 
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This  windowing  is  useful  if  only  a  part  of  the  seismogram  corresponding  with 
certain  ray  phase  velocities  is  required  is  desired.  When  this  is  used,  the 
upper  limit  in  the  integrand  is  defined  by  c^,,,  and  the  parameter  FAC  and 
the  asymptotic  integral  technique  is  not  used.  The  use  of  the  Hankel,  rather 
than  Bessel,  function  is  recommended  when  phase  velocity  filtering  is  per¬ 
formed.  The  filtering  will  cause  numerical  noise,  but  will  also  have  the  advan¬ 
tage  that  the  upper  limit  of  wavenumber  integration  will  be  reduced,  and 
hence  the  computations  will  be  accomplised  faster. 

Examples 

To  illustrate  these  concepts,  consider  the  RDD  Green’s  function  for  a 
wholespace,  (lb)  and  (3b).  The  model  consists  of  P-velocity  of  6.0  km/s,  an  S- 
velocity  of  3.5  km/s,  and  a  density  of  2.8  gm/cm3.  The  source  is  5  km  beneath 
the  receiver,  and  FAC  =  3.0.  A  128  point  time  series  is  generated  with  a  sam¬ 
pling  interval  of  1.0  sec,  and  a  triangular  pulse  with  r  =  2.0  sec  is  used  to  gen¬ 
erate  velocity  time  histories  in  units  of  cm/s  for  a  source  moment  of  1020  dyne- 
cm,  at  distances  of  1  -  500  km.  The  time  domain  damping  parameter  a  -  0. 02 
so  that  the  periodicity  effect  is  reduced  by  a  factor  of  exp(~2. 54).  The  RDD 
Green’s  function  is  a  good  choice  to  illustrate  the  problems  of  numerical  noise 
because  the  and  S-wave  arrival  falls  off  more  rapidly  with  distance  than  the 
P-wave  arrival. 
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Figure  4  shows  the  time  histories  when  L  =  250km  is  used.  The  Bou- 
chon  (1981)  box  bounded  in  distance  by  L/2  and  in  time  by  L/2Vp  is  shown 
by  the  long  dashed  lines.  This  rectangular  box  is  used  since  one  often  com¬ 
putes  a  record  section  rather  than  a  single  trace,  and  thus  the  size  of  the 
acceptable  time  window  is  constrained  by  the  distance  trace.  The  direct  and 
wavenumber  integration  periodicity  arrivals  are  indicated  by  the  light  solid 
and  short  dashed  lines,  with  the  latter  indicating  those  arrivals  that  have 
been  wrapped  around  due  to  the  temporal  periodicity  properties.  The  extent 
of  the  periodicity  problems  is  readily  apparent,  as  is  the  fact  that  the  signal 
within  the  Bouchon  box  is  relatively  clean.  Part  of  the  "simplicity"  of  the 
waveforms  is  due  to  the  fact  that  the  choice  of  a  reduced  the  amplitude  of  the 
arrivals  folded  in  time. 

Figure  5  shows  the  same  model,  except  that  L= 500km.  Because  of  the 
change  in  L,  the  Bouchon  box  becomes  larger,  and  the  seismograms  are  given 
to  a  larger  distance.  A  P-wave  noise  arrival  due  to  the  temporal  and  spatial 
periodicity  appears  prior  to  the  direct  P-wave  arrival.  Figure  5  also  shows 
some  of  the  k  =  0  noise  at  distances  greater  than  300  km  at  an  arrival  time  of 
about  1  sec.  This  noise  (Herrmann  and  Mandal,  1986)  was  reduced  using  the 
modified  integration  rule,  but  is  a  low  frequency  arrival,  and  is  reduced  fur¬ 
ther  by  making  L  smaller.  Figure  6  is  a  similar  plot,  but  uses  a  reduction 
velocity  of  6.0  km/sec.  The  reduced  travel  time  will  provide  a  cleaner  begin¬ 
ning  of  the  trace  history,  but  will  not  give  any  additional  respite  from  the 
Bouchon  condition. 

Figure  7  shows  the  effect  of  using  the  Hankel  function  rather  than  the 
Bessel  function  together  with  a  phase  velocity  filter.  The  phase  velocity  filter 
window  of  [  16,  8,  2,  1  ]  km/s  was  appropriate,  since  at  large  offset,  the  rays 
corresponding  to  P-  and  S-wave  arrivals  would  propagate  at  6.0  amd  3.5  km/s 
respectively.  The  Hankel  function  was  used  in  the  integration  when  its  argu¬ 
ment,  e.g.,  kr,  was  greater  than  6.0.  Otherwise  the  Bessel  function  is  used. 
Since  this  can  generate  some  low  frequency  noise  near  k  **  0,  the  time  domain 
damping  parameter  was  set  to  a  =  0.01  rather  than  the  or  =  0.02  used  in  the 
preceding  figures.  This  reduces  the  amplitude  of  late  arriving  noise,  but  also 
increases  the  amplitude  of  the  time  domain  periodicity  noise. 

There  are  two  impressive  points  that  can  be  made  comparing  Figures  6 
and  7.  First,  the  wavenumber  integration  synthetics,  left  side  of  Figure  7, 
agree  very  well  with  the  analytic  time  histories,  Figure  7  right,  in  the  dis¬ 
tance  range  of  25  -  350  km,  and  in  the  time  window  up  to  83  seconds  after  the 
P-wave.  This  is  significantly  better  in  time  and  distance  than  the  Bouchon 
window  shown  in  Figure  6.  This  is  in  spite  of  the  fact  that  the  computational 
effort  is  roughly  the  same  since  in  this  case  FACkVniiB  ~  <a/ 1.  The  implication 
is  that  if  the  range  of  interest  were  1  -  250  km  and  the  time  window  were  still 
41.67  sec,  e.g., the  Bouchon  box,  then  it  should  be  possible  to  use  a  smaller 
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Fig.  4.  HDD  Green’s  functions  for  L  =  250  km.  All  traces  start  at  -2.0  sec  and  end  at  125  sec. 
The  distances  and  peak  velocities  (cm/s)  are  given  to  the  right  of  each  trace.  The  long  dashed 
box  shows  the  Bouchon  window,  the  light  solid  lines  indicate  the  direct  arrivals  and  those 
expected  from  periodicity  in  space  due  to  wavenumber  integration,  and  the  light  short- 
dashed  lines  the  folding  of  these  arrivals  due  to  the  temporal  periodicity  due  to  the  use  of  a 
discrete  Fourier  transform. 


value  of  L  to  not  only  increase  the  time  window  but  to  also  reduce  the  eompu 
tation  time. 
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Fig.  5.  RDD  Green’s  functions  for  L  =  500  km.  All  traces  start  at  -2.0  sec  and  end  at  125  sec. 
The  distances  and  peak  velocities  (cm/s)  are  given  to  the  right  of  each  trace.  The  long  dashed 
box  shows  the  Bouchon  window,  the  light  solid  lines  indicate  the  direct  arrivals  and  those 
expected  from  periodicity  in  space  due  to  wavenumber  integration,  and  the  light  short- 
dashed  lines  the  folding  of  these  arrivals  due  to  the  temporal  periodicity  due  to  the  use  of  a 
discrete  Fourier  transform. 


Guidance  in  Choosing  Correct  Parameters 

If  high  frequency  synthetic  seismograms  are  desired,  computer  runs  will 
necessarily  be  lengthy.  The  criteria  for  choosing  the  time-domain  damping 
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Fig.  6.  RDD  Green’s  functions  for  L  =  500  km.  A  reduced  travel  time  plot  is  used.  The  start 
time  of  each  trace  is  r/6.0  -  2.0  seconds  after  the  origin  time;  r  is  the  epicentral  distance  in 
km.  All  traces  are  123  seconds  long.  All  traces  start  at  The  distances  and  peak  velocities 
(cm/s)  are  given  to  the  right  of  each  trace.  The  long  dashed  box  shows  the  Bouchon  window, 
the  light  solid  lines  indicate  the  direct  arrivals  and  those  expected  from  periodicity  in  space 
due  to  wavenumber  integration,  and  the  light  short-dashed  lines  the  folding  of  these  arrivals 
due  to  the  temporal  periodicity  due  to  the  use  of  a  discrete  Fourier  transform. 


parameter  a  is  straightforward.  The  choice  of  the  wavenumber  sampling 
parameter  L  is  not  as  simple.  One  useful  technique  used  to  settle  upon  a 
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Fig.  7.  Comparison  of  wavenumber  integration  (left)  and  analytic  (right)  RDD  Green’s  func¬ 
tions.  For  the  wavenumber  integration  synthetics,  the  phase  velocity  filter  was  defined  by  the 
window  [  16,  8,  2,  1],  the  Hankel  function  was  used,  and  L  =  500  km.  For  both  displays,  a 
reduced  travel  time  plot  is  used.  The  start  time  of  each  trace  is  r(6.0  -  2  0  seconds  after  the 
origin  time;  r  is  the  epicentral  distance  in  km.  All  traces  are  123  seconds  long.  All  traces 
start  at  The  distances  and  peak  velocities  (cm/s)  are  given  to  the  right  of  each  trace.  In  the 
wavenumber  integration  synthetics,  the  light  solid  lines  indicate  the  direct  arrivals  and 
those  expected  from  periodicity  in  space  due  to  wavenumber  integration,  and  the  light  short- 
dashed  lines  the  folding  of  these  arrivals  due  to  the  temporal  periodicity  due  to  the  use  of  a 
discrete  Fourier  transform.  Since  the  Hankel  function  is  used,  no  wavenumber  integration 
noise  propagates  with  negative  velocity. 


suitable  value  is  to  keep  a  and  NAt  fixed,  whil°  making  N  small.  For  the  Bou- 
chon  sampling,  decreasing  N  by  a  factor  of  2,  should  decrease  the  number  of 
samples  in  the  k  -  to  plane  by  a  factor  of  4.  In  addition,  the  resultant  time  his¬ 
tories,  consisting  of  low  frequency  components,  will  exhibit  the  same  sensitiv¬ 
ity  to  the  choice  of  L  as  the  desired  high  frequency  time  histories.  Thus  if  the 
low  frequency  time  history  looks  good  in  terms  of  the  wavenumber  integra¬ 
tion  noise,  the  high  frequency  time  histories  also  will  look  good. 

NUMERICAL  ISSUES  IN  EVALUATING  PROPAGATOR  MATRICES 

The  integrands  in  (1)  are  highly  dependent  upon  the  difference  between 
the  source  and  receiver  depths,  and  also,  in  the  case  of  a  layered,  medium, 
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the  layer  thicknesses.  For  a  layered  medium  the  propagator  matrices  contain 
hyperbolic  sine  and  cosine  terms  when  k>kv,  where  V  is  either  the  P-  or  S- 
wave  velocity.  At  high  frequencies  or  large  values  of  vvd,  where  d  is  the  layer 
thickness,  the  size  of  the  hyperbolic  terms  can  quickly  exceed  the  size  of  float¬ 
ing  point  numbers  in  a  computer. 

This  problem  is  handled  very  effectively  by  introducing  an  extended 
floating  point  notation.  The  SH  propagator  matrix,  given  in  the  Appendix,  is 
of  the  form: 


A(z)  = 


cosh  X  b  sinh  X 
asinhX  coshX 


For  large  values  of  X,  this  can  be  rewritten  as 


A(z)  -  g 


1  +  e 


-2X 


b  1  -  e 


-2X 


a  1-e 


-2X 


1  +  e 


2X 


Note  that  now  the  exponential  terms  within  the  matrix  will  now  only  under¬ 
flow,  and  can  be  set  to  zero  when  X  is  sufficiently  large.  However,  the  leading 
ex  term  is  not  multiplied,  but  the  X  value  is  saved.  The  A(z)  matrix  will  now 
look  like 


A(z)  =  exA'(z) 


and  matrix  multiplication  of  two  such  propagators  will  give 

A(zi)A(z2)  =  eXl+X2A'(z1)A'(z2) 


Computationally  a  "modified  matrix  product"  is  saved  as  are  the  £Xj. 
Since  the  expression  for  receiver  displacements  for  SH  waves  (38)  is  a  ratio  of 
two  sets  of  propagator  matrices,  the  resultant  expression  will  look  like 

U,  =  e^NUMi~^'DEN,well_defined_number. 

If  the  exponential  term  is  too  small,  U,  is  set  to  zero. 

This  artifice  works  well.  However,  there  is  an  assumption  inherent  in 
the  use  of  (38)  or  (19)  or  (20),  and  that  is  that  the  exponential  eigenfunctions 
always  decrease  from  the  source  to  the  receiver,  since  the  propagators  are 
taken  from  the  base  of  the  layer  stack  to  the  surface  and  the  above  factoriza¬ 
tion  assumes  that  they  increase  upwards.  Thus  there  may  be  some  problems 
with  sme  models  having  low  velocity  zones. 

Another  computational  problem  arises  with  the  use  of  propagator  matri¬ 
ces.  Harkrider  (personal  communication)  mentioned  the  inherent  instability 
of  the  compound  P-SV  matrices  at  low  frequencies.  While  this  may  be  a 
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problem  in  single  precision,  we  have  not  found  it  to  be  so  when  computations 

are  performed  in  double  precision. 

REFERENCES 

Apsel,  R.  J.,  and  J.  E.  Luco  (1983).  On  the  Green’s  functions  for  a  layered 
half-space.  Part  II.,  Bull.  Seism.  Soc.  Am.  73, 931-951. 

Bouchon,  M.  (1981).  A  simple  method  to  calculate  Green’s  functions  for  elas¬ 
tic  layered  media,  Bull.  Seism.  Soc.  Am.  71, 959-971. 

Brigham,  E.  O.  (1974).  The  Fast  Fourier  transform,  Prentice- Hall,  Engle¬ 
wood  Cliffs. 

Fuchs,  K.,  and  G.  Muller  (1971).  Computation  of  synthetics  seismograms 
with  the  reflectivity  method  and  comparison  with  observations,  Geophys. 
J.  23,417-433. 

Haskell,  N.  A.  (1963).  Radiation  patters  of  Rayleigh  waves  from  a  fault  of 
arbitrary  dip  and  direction  of  motion  in  a  homegeneous  medium,  Bull. 
Seism.  Soc.  Am.  53,  619-642. 

Haskell,  N.  A.  (1964).  Radiation  pattern  of  surface  waves  from  point  sources 
in  a  multi-layered  medium,  Bull.  Seism.  Soc.  Am.  54,  377-393. 

Herrmann,  R.  B.  (1975).  A  student’s  guide  to  the  use  of  P  and  S  wave  data 
for  focal  mechanism  determination,  Earthquake  Notes  46,  29-40. 

Herrmann,  R.  B.,  and  B.  Mandal  (1986).  A  study  of  wavenumber  integration 
techniques,  Earthquake  Notes  57, 33-40. 

Herrmann,  R.B.,  and  C.  Y.  Wang  (1985).  A  comparison  of  synthetic  seismo¬ 
grams,  Bull.  Seism.  Soc.  Am.  75,  41-56. 

Hudson,  J.  A.  (1969).  A  quantitative  evaluation  of  seismic  signals  at  teleseis- 
mic  distances  -  I.  Radiation  from  point  sources.  Geophys.  J.  Roy.  astr. 
Soc.  18,233-239. 

Kennett,  B.  L.  N.  (1983).  Seismic  Wave  Propagation  in  Stratified  Media, 
Cambridge  University  Press,  Cambridge. 

Mallick,  S.,  and  L.  N.  Frazer  (1987).  Practical  aspects  of  reflectivity  model¬ 
ing,  Geophysics  52,  1355-1364. 


-46- 


Phinney,  R.  A.  (1965).  Theoretical  calculation  of  the  spectrum  of  first  arrivals 
in  layered  elastic  mediums,  J.  Geophys.  Res.  70,  5107-5123. 

Pujol,  J.,  and  R.  B.  Herrmann  (1990).  A  student’s  guide  to  point  sources  in 
inhomogeneous  media,  Seism.  Res.  Letters  61,  209-224. 

Saikia,  C.  K.  (1993).  Modified  frequency-wavenumber  algorithm  for  regional 
seismograms  using  Filon’s  quadrature  -  Modeling  of  Lg  waves  in  eastern 
North  America,  Geophys.  J.  Int.  (in  review). 

Wang,  C.  Y.,  and  R.  B.  Herrmann  (1980).  A  numerical  study  of  P-,  SV-,  and 
SH-wave  generation  in  a  plane  layered  medium.  Bull.  Seism.  Soc.  Am. 
70, 1015-1036. 

Watson,  T.  H.  (1970).  A  note  on  fast  computation  of  Rayleigh  wave  dispersion 
in  the  multi-layered  elastic  half-space,  Bull.  Seism.  Soc.  Am.  60,  161- 1 66. 


-47- 


APPENDIX 


P-SV  Propagator 

The  components  of  the  4x4  P-SV  propagator  from  Haskell  (1964)  are  as 
follow: 


au  =  y  cosh  vaz  ~(y  -  l)cosh  vez 

ai 2=  -(y-l)sinhvaz/va  +  yv^sinh^z/k2 

a13  =  -  (cosh  vaz- cosh  v0z)f p 


aX4  = 


(k2 sinh  vaz!va-ve  sinh  vfiz ) j/ p 


a2i  =  yvasinh  vaz-k2(y- l)sinh  v0z!vp 
a22=  -  (y  - 1)  cosh  vuz  +  y  cosh  vpz 


a2 3=  -  va  sinh  v^z  +  k2  sinh  vpz/vp\/p 


a24  -  ~k2a13 

a3i  =  py(y~  l)(cosh  v^z-cosh  vpz) 


a32  =  p[^~  (y  -  l)2  sinh  vailva  +  y2  v0  sinh  vpz  / k2  j 


a33  ~  a22 

a34  =  ~  k2a12 

a4i  =  PY2va  s*nh  v0z/k2  -  p(y-l)2  sinh  vpz!vp 

a42  =  ~  a3i  ^k2 
a43  =  a2 1  /  k2 
a44  =  3 1 1 

The  compound  matrix  of  a  4x4  matrix,  a^  is  a  6x6.  However,  the  struc¬ 
ture  of  the  compound  matrix  A I  ^,  =  aikaji  -  a^a^,  permits  a  reduction  to  a  5x5 
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matrix,  (Watson,  1970),  which  reduces  the  number  of  matrix  multiplications 
required.  Associating  compound  matrix  doublets  { 12,  13,  14,  23,  24,  34  ), 
with  the  indices  {1,  2,  3,  4,  5,  6  }  ,  and  defining  Ay  =  Al[]>,  using  the  above 
mapping,  we  see  that 

^3j  =  -  A»j  /  k2  ,  j  *  3, 4 

Ai3  =  -k2Aj4  ,  j  *■  3, 4 


A33  -  A 44  , 

A34=-(A44-l)/k2 

In  addition,  Ay  =  A7^  7_j  for  i  =  2  to  6  and  j  =  8  -  i  to  6. 

When  this  matrix  is  multiplied  by  a  vector,  e.g.,  the  elements  of  the 
compound  H  matrices,  we  note  that  the  third  element  of  the  product,  G3,  is 
related  to  the  fourth  element,  G4 ,  by  the  relation  G3  =  -  k2G4 .  For  the  com¬ 
pound  H  matrix,  the  third  and  fourth  elements  are  related  by  H4  =-k2H3. 
The  1x6  and  6x6  multiplication,  GTA  yields  the  same  results  as  multiplying  a 
1x5  by  a  5x5,  G'  A'  if  we  do  the  following  to  the  G  vector  and  the  A  matrices: 

a)  Drop  the  third  element  of  the  compount  G. 

b)  Drop  the  third  row  and  column  of  A  to  form  the  initial  5x5. 

c)  Multiply  the  off  diagonal  elements  of  the  third  row  of  the  new  A'  by  2. 

d)  Define  the  third  diagonal  element  A33  =  2A44  -  1. 

Thus  the  original  6x6  compound  matrix 


an 

a12 

a13 

a14 

a  15 

a  16 

a21 

a22 

a23 

a24 

a25 

a26 

a31 

a32 

a33 

a34 

a35 

a36 

a41 

aj2 

a43 

a44 

a45 

a46 

a51 

a52 

a63 

a64 

a55 

a  56 

.^l 

362 

363 

3fi4 

365 

366 

becomes 


3U 

a12 

a  14 

a15 

a16 

a21 

a22 

a24 

a25 

3  26 

2341 

2342 

2a44  - 1 

2a  45 

2a46 

a61 

a52 

a  54 

a65 

3  56 

3fil 

3fi2 

a  64 

a65 

3fi6 

The  original  compound 


-49- 


G  =|gn.gi2.gl3,gu.gl6.gl6j 

becomes 

G'  =  j^gll>  gl2»  gl4>  gl5,gl6^J  1 

and  the  original  compound 

HT  =  ^h11,h2i,h31,h41th6i,h61j 

becomes 

H'7  =^h11,h2i,2h4i,h6iih61j  . 

We  will  compute  the  propagators  using  the  primed  matrices.  Tb  evaluate  (19) 
or  (20),  we  can  quickly  define  the  needed  unprimed  vector  from 


>f 

‘  ri  ' 

r2 

*2 

r3 

-k2r^ 

r4 

r3 

r5 

ri 

-re. 

_  rs  . 

The  elements  of  the  modified  5x5  compound  propagator  matrix  are  as 
follow: 

An  =  CPCQ  -  (A33  - 1)/2 
A'i2  =  ( -  CQX  +  k2CPY)  /  p 

A'13  =  -  [(2y  -  IX 1  -  CPCQ)  +  yXZ /  k2  +  (y  -  l)k2WY] /  p 

A'14  =  (CPZ-k2CQW)/p 

A'15  =-[2(1-  CPCQ)k2  +  WYk4  +  XZ]  /  p2 

Ak  =  p[-(r~  l)2CQW+y2CPZ/k2] 

A22  =  CPCQ 

A£3  =  -[(y-l)CQW-yCPZ/k2] 

A^,=  -WZ 
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A25  -  A'14 

A31  =  2 p[y(y -  1X2 y -  1X1  - CPCQ)  +  r3XZ/k2  +  (y~  l)3k2WY] 
A^2  =  2[rCQX-k2(r-l)CPY] 

A'33  i  +  2[2y(y  -  1X1  -  CPCQ)  +  y  2XZ  /  k2  +  (y  -  l)2k2WY] 

A34  —  -  2k2  A23 
A35  «  —  2k2A'13 

Aii  =  p[-y2CQX/k2  +  (y-  1)2CPY) 

a;2  =  -xy 

A43  =  “  2  A32/k2 
A44  =  A22 
A46  ~  Aj2 


Abi  ~  —  p‘ 


(y  -  m  2(  1  -  CPCQ)y2  +  (y - 1)2 WYk2 1+  y4X Z/k2 


/k2 


Am  -  Aii 
A'63=-|Ai1/k2 


where  we  define 


CPCQ  =  COSP  *  COSQ 

CPY  =  COSP  *  Y  XY  =  X  *  Y 

CPZ  =  COSP  *  Z  XZ  =  X*Z 

CQW  =  COSQ  *  W  WY  =  W  *  Y 
CQX  =  COSQ  *  X  WZ  =  W  *  Z 
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COSP  =  cosh  vaz  X  =  va  sinh  vaz 
COSQ  =  coshv^z  Y  =  sinh^z/^ 
W  =  sinh  vaz /  va  Z  =  sinh  vfiz 


For  a  fluid  layer,  the  values  of  the  vertical  displacement  and  stress  must 
be  propagated,  but  reduce  the  need  for  complicated  if  constructs  in  the  pro¬ 
gram,  the  radial  stress  and  displacement  values  are  propagated  unchanged 
in  this  modified  matrix.  The  true  2x2  propagator  matrix  for  a  fluid  consists  of 
the  inner  2x2  submatrix.  The  propagator  used  is  as  follows: 


a(z)  = 


0 

cosh  vaz 
-psinh  vaz  /  va 
0 


0  O' 

-va  sinhvaz/p  0 
cosh  vaz  0 

0  1 


The  corresponding  5x5  modified  compound  propagator  matrix  is 


cosh  vaz 
-  p  sinh  vaz/v0 
A(z)  =  0 

0 


-va  sinhvaz/p  0 

cosh  vaz  0 

0  1 

0  0 


cosh  vaz 


-psinh  vazlva 


0  -  va  sinh  vazt  p  coshvaz 


SH  Propagator 

The  components  of  the  2x2  SH  propagator  follow  follow  the  Herrmann 
(1979)  modification  of  the  Haskell  (1964)  matrices  to  factor  out  an  appraent 
singularity  at  ca  ~  0.: 


A(z)  = 


coshv^z  sinh  vfiz/ p(i2Vp 
pfPvp  sinh  vpz  cosh  vfiz 


If  the  top  or  bottom  layers  of  the  structure  are  fluid,  a  pseudo- 
propagator  is  introduced  for  SH  to  avoid  the  use  of  complicated  if  structures 
in  the  source  program  with  the  following  matrix: 


A(z)  = 


Functional  Form  of  Integrands  as  <v  ->0ork»  kVm.n 

Explicit  expressions  for  the  earthquake  and  explosion  double  couples  for 
surface  receiver  and  a  buried  source  in  a  halfspace  were  given  by  Herrmann 
and  Wang  (1985).  This  required  care  in  taking  the  limits  of  all  the  terms  in 
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the  specific  expressions. 

For  each  Green’s  functions,  we  may  wish  to  evaluate  an  integral  of  the 

form: 


JnKm  =  J  kme_khJn{kr)dk  . 

o 

Specific  forms  can  be  found  by  integration  and  differentiaion  of  the  Sommer- 
feld  integral  with  respect  to  r  or  z. 

For  the  Green’s  functions,  the  asymptotic  form  was  evaluated,  and  since 
only  two  wavenumber  values,  kj  and  k2,  are  used,  two  constants  are  defined 
needed.  In  addition,  F^  and  F16  vary  as  k'1  for  large  k.  To  handle  this  case 
stably  the  amsyptotic  fit  is  made  to  kF12  and  to  kF16.  Table  A.l  shows  the 
particular  asymptotic  coefficients  evaluated.  The  lowest  order  term  is  k  is 
that  which  results  from  taking  the  limits.  By  combining  this  table  with  the 
expresion  in  (1),  the  necessary  integrals  are  defined. 
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Table  A.  1.  Functional  fit  to  asymptotic  trend 
Fj  Asymptotic  Fit  Function  Fit 

Fx  (Bk+Ck2)e-kh  Fx 

F2  (A+Bkkf"1  F2 

F3  (Bk+CkV1*  F3 

F4  (A+Bkle^  F4 

F6  (Bk+Ck2)e_kh  Fs 


f6 

(A+Bkle^ 

f6 

f7 

(Bk+Ck2)e_kh 

f7 

f8 

(A+Bkle-1* 

f8 

f9 

(Bk  +  Ck2)e-kh 

f9 

F10 

(A+Bk)e"*h 

Fio 

Fu 

(Bk  +  Ck2)e-kh 

Fn 

f12 

(A+Bk)e_kh 

kF12 

f13 

(A+Bkle^ 

f13 

f14 

(A+Bk)e_kh 

f14 

f15 

(A+Bkle"1* 

kF16 

f16 

(Bk  +  Ck2)e"kh 

f16 
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The  necessary  integrals  are 

n  m  JnKm 


0  0 
0  1 
0  2 
0  3 

1  -1 
1  0 
1  1 
1  2 

1  3 

2  0 

2  1 
2  2 
2  3 


* 

R3 


R® 

6z3  -  9ZT2 
R7 


R2  =  r2  +  h2 
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SPECTRAL  EXAMINATION  OF  THE  16  JUNE  1992  EARTHQUAKE 
AND  QUARRY  BLAST  NEAR  EVANSVILLE,  INDIANA 

K.D.  Hutchenson  and  R.B.  Herrmann 

ABSTRACT 


On  16  June  1992,  an  2.3  earthquake  occurred  in  southwest¬ 
ern  Indiana,  near  Evansville.  This  area  is  part  of  the  Dlinois 
Basin  coal  belt,  an  area  of  active  surface  mines  with  numerous 
strip-mine  blasts  daily.  The  co-location  of  earthquakes  and 
strip-mine  blasts  enable  spectral  comparisons  without  signifi¬ 
cant  concern  for  differences  due  to  path  propagation  effects. 

Discriminating  between  the  two  types  of  events  can  be  done 
visually  due  to  the  distinctive  appearance  of  the  Rg  phase  in 
strip-mine  blasts  and  high  frequency  coda  of  earthquakes.  A 
strong  Rg  phase  is  indicative  of  shallow  source  depths.  How¬ 
ever,  earthquakes  previously  located  at  shallow  depths  else¬ 
where  within  the  Illinois  Basin  do  not  exhibit  a  distinctive  Rg 
phase,  indicating  either  poor  control  in  focal  depth  determina¬ 
tion  or  a  fundamental  difference  in  source  mechanism.  Visual 
and  spectral  examination  shows  that  earthquakes  are  richer  in 
energy  at  higher  frequencies  than  strip-mine  blasts.  Earth¬ 
quakes  have  significant  energy  at  20-30  Hz,  while  the  significant 
energy  content  of  blasts  is  closer  to  10  Hz.  The  significant  differ¬ 
ence  compared  to  previous  earthquake-nuclear  explosion  dis¬ 
criminant  studies  is  that  the  chemical  explosion  has  reduced 
high  frequency  content  compared  to  the  earthquake. 

INTRODUCTION 

On  16  June  1992,  an  m^  2.3  earthquake  occurred  in  southwestern  Indi¬ 
ana  near  Evansville.  The  earthquake  was  recorded  by  the  Central  Missis¬ 
sippi  Valley  Seismic  Network  (CMVSN)  (Figure  1).  Although  this  area  is  not 
as  seismically  active  as  the  New  Madrid  Seismic  Zone  to  the  southwest,  small 
earthquakes  (1.7  <  m^  <  2.9)  are  not  infrequent,  averaging  one  instrumen- 
tally  located  event  every  thirteen  months  since  September  1982.  This  study 
was  prompted  by  the  fortuitous  co-location  of  the  16  June  earthquake  (only 
the  third  digitally  recorded  earthquake  located  in  this  area)  with  a  nearby 
strip-mine  blast  of  similar  size  (m^  2.65)  which  occurred  just  hours  after  the 
earthquake  (Figure  2). 


Fig.  1.  Seismicity  map  for  the  New  Madrid  Seismic  Zone,  Illinois  Basin,  and  surrounding 
areas  showing  the  CMVSN  stations.  Two  insets  show  the  Illinois  Basin  outline  (after 
Collinson  et  ai,  1988)  and  the  area  of  study  in  eastern  Illinois  and  southwestern  Indiana  (cir¬ 
cles  (o)  show  prior  seismicity,  with  the  two  squares  (□)  denoting  the  1  May  1985  and  16  June 
1992  events;  asterisks  (♦)  show  two  strip  mine  blasts  at  the  same  location. 


The  primary  purpose  of  this  paper  is  to  examine  the  distinguishing  fea¬ 
tures  between  the  earthquake  (19:28:49.0  UT)  and  strip-mine  blast 
(23:11:13.6  UT)  of  16  June  1992  (Figure  1)  on  the  basis  of  several  attributes, 
including  spectral  differences.  In  addition,  these  results  will  be  compared 
with  similar  events,  both  in  the  immediate  area  and  elsewhere  within  the 
Illinois  Basin.  It  will  be  shown  that  differentiating  between  these  earth¬ 
quakes  and  strip-mine  blasts  within  the  Illinois  Basin  is  rather  trivial. 

The  intent  of  this  study  is  not  to  blindly  apply  regional  nuclear  explo¬ 
sion-earthquake  discrimination  results  from  other  areas  {e.g.  Pomeroy  et  al., 
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Fig.  2.  WDIN  analog  record  segment  of  16  June  1992  showing  the  earthquake  and  several 
strip-mine  blasts.  The  events  with  the  low  frequency  ’tail"  (coda)  are  strip-mine  blasts  (three 
within  50  km;  three  more  farther  away).  As  shown  by  this  seismogTam  segment,  strip  mine 
blasts  and  earthquakes  are  quite  different  in  appearance,  primarily  in  the  appearance  of  the 
surface  wave  with  respect  to  the  body  wave  arrivals.  The  majority  of  the  "events"  on  this 
record  are  associated  with  vehicular  traffic. 


1982;  Bennett  and  Murphy,  1986;  Taylor  et  al.,  1989),  but  to  examine 
observed  characteristics  which  may  aid  in  separating  these  two  events.  The 
differences  are  due  to  the  nature  of  the  source;  low  yield  nuclear  explosions 
are  considered  point  sources  while  large  strip-mine  blasts  are  a  series  of 
sources.  However,  this  area  would  provide  a  good  test  case  for  more  detailed 
work,  in  that  events  have  both  similar  size  (m1/K)  and  location,  except  for 
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depth.  Discrimination  studies  undertaken  elsewhere  within  the  Dlinois 
Basin  and  the  New  Madrid  seismic  zone  have  examined  the  classical  discrim¬ 
inants  (Amjad,  1991). 

DESCRIPTION  OF  THE  STUDY  AREA 

The  Illinois  Basin  (600  km  by  320  km)  is  located  primarily  in  Illinois, 
extending  into  Indiana,  and  Kentucky  (Collinson  et  ai,  1988)  (Figure  1).  The 
Paleozoic  sedimentary  fill  (primarily  carbonates  with  lesser  amounts  of  silt- 
stone/shale  and  sandstones)  range  in  age  from  Early/Middle  Cambrian  to 
Early  Permian.  Parts  of  the  basin  are  capped  by  Pleistocene  loess  deposits  of 
varying  thickness  (Willman  and  Frye,  1970).  The  thickest  sediments  are 
along  the  southern  border  of  the  basin,  in  central  and  southern  Illinois, 
southwestern  Indiana,  and  western  Kentucky,  thinning  both  depositionally 
and  erosionally  towards  the  basin  boundaries. 

The  Wabash  River  (boundary  between  Illinois  and  Indiana)  area  has 
been  examined  in  prior  geological  and  geophysical  studies  (Bristal  and 
Buschbach,  1971;  Bristal  and  Treworgy,  1979;  Ault  and  Sullivan,  1982;  Ault 
et  ai,  1985;  Nelson  and  Bauer,  1987).  Many  high-angle  normal  faults  have 
been  mapped  trending  NNE,  subparallel  to  the  Wabash  River.  Offsets  along 
the  faults  may  diminish  with  depth  (Sexton  et  al.,  1988). 

Seismic  activity  has  been  reported  in  the  eastern  Illinois,  western  Indi¬ 
ana  area  since  the  1800’s  (NuttU,  1983).  The  Wabash  River  Seismic  Zone 
(WRSZ)  was  identified  as  a  result  of  the  seismic  activity  and  the  many 
mapped  faults  (Nuttli  and  Brill,  1981).  The  WRSZ  has  been  considered  a 
source  zone  in  terms  of  hazard  analysis  (Barstow  et  al.,  1981).  However, 
recent  work  (Taylor,  1991)  has  argued  that  the  WRSZ  is  not  a  separate  zone 
and  should  be  included  together  with  eastern  Illinois  as  one  zone.  Seismic 
activity  in  eastern  Illinois,  including  the  largest  event  recorded  by  the 
CMVSN  (Taylor  et  ai,  1989),  may  be  associated  with  basement  structures  in 
the  region  (Hamburger  and  Rupp,  1988;  Pratt  et  ai,  1989;  Taylor,  1991). 

In  addition  to  the  historical  and  instrumentally  located  seismic  activity 
along  the  Wabash  River,  there  are  numerous  strip  mine  blasts.  Surface 
mines  are  quarrying  coal  from  Pennsylvanian  deposits  (Hasenmueller  and 
Carr,  1983;  Hasenmueller  and  Wiegand,  1980  (rev.  1985);  Harper,  1985). 
Many  artificial  events  are  recorded  by  the  CMVSN  on  a  weekly  basis. 

DATA  PROCESSING  AND  RESULTS 

Figure  1  shows  the  seismic  activity  in  the  New  Madrid  Seismic  Zone 
and  Illinois  Basin  area.  Insets  show  the  study  area  in  southwestern  Indiana 
and  eastern  Illinois,  with  the  located  seismicity  and  several  strip-mine  blasts. 
Also  shown  are  the  CMVSN  stations,  especially  WDIN  and  BPIL,  which  are 
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used  in  the  following  analysis.  WDIN  (Wadesville,  IN),  is  located  approxi¬ 
mately  40  km  from  the  active  strip  mines,  while  BPIL  (Belle  Praire,  IL)  is 
farther  away,  nearly  100  km  from  the  active  mines. 


Table  1 

Station  Locations 


Station 

Location 

Lat 

(°N) 

Lon 

(°W) 

BPIL 

Belle  Praire,  IL 

38.20 

88.59 

CCMO 

Creve  Coeur,  MO 

38.72 

90.47 

CSIL 

Creal  Springs,  IL 

37.63 

88.79 

CIRL 

Cave  In  Rock,  IL 

37.51 

88.11 

DON 

Dongola,  MO 

37.18 

89.93 

ELC 

Elco,  IL 

37.28 

89.23 

FVM 

FYench  Village,  MO 

37.98 

90.43 

GOIL 

Rosebud,  IL 

37.29 

88.58 

NHIL 

New  Haven,  IL 

37.93 

88.17 

SLM 

St.  Louis,  MO 

38.64 

90.24 

TYS 

Tyson  Valley,  MO 

38.53 

90.57 

WDIN 

Wadesville,  IN 

38.09 

87.72 

WDIL 

West  Salem,  IL 

38.50 

88.07 

While  several  permanent  CMVSN  stations  are  located  along  the 
Wabash  River,  most  lie  on  the  Illinois  side  (Table  1).  The  network  has  been  in 
operation  since  1974.  All  sites  have  vertical  component  seismometers 
(L4-C’s)  except  for  FVM  and  SLM  which  have  three  component  Benioff  seis¬ 
mometers.  Data  are  telemetered  to  St.  Louis  University  where  they  are  digi¬ 
tized  at  100  Hz. 

Table  2  lists  the  parameters  of  all  seismicity  recorded  in  the  study  area, 
in  addition  to  several  strip-mine  blasts.  Only  three  of  the  earthquakes  were 
digitally  recorded.  Of  these,  two  are  used  in  this  study;  the  13  February  1987 
event  is  virtually  unusable  due  to  poor  signal-to-noise. 

Two  digitally  recorded  earthquakes  (1  May  1985  and  16  June  1992)  and 
two  recent  strip-mine  blasts  of  similar  size  are  examined  in  this  study.  Each 
event  was  located  using  a  version  of  FASTHYPO  (Herrmann,  1979a)  with  the 
UPLANDS  velocity  model  (Stauder  et  al.,  1991)  and  available  CMVSN  sta¬ 
tions  (Table  1  and  Figure  1).  In  addition,  station  BLO  (Bloomington,  IN) 
reported  phase  data  for  several  of  the  events  (M.  Hamburger,  pers.  comm.),  as 
did  several  Kentucky  stations  (LLKY,  MOKY,  SMKY,  and  SOKY)  (R.  Street, 
pers.  comm.),  aiding  in  their  location,  especially  with  depth.  Due  to  the  sta¬ 
tion  geometry,  error  ellipses  for  all  events,  both  blasts  and  earthquakes  in  the 
Wabash  River  Valley,  are  generally  oriented  in  a  NE-SW  direction. 
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Table  2 


Earthquakes  and  strip  mine  blasts  in  southwestern  Indiana 


Date 

Time) 

(UT) 

Lat 

(°N) 

Lon 

(°W) 

Depth 

(km) 

mb 

comments 

2  Sep  1982 

11:30:59.3 

37.95 

87.75 

12.0 

1.3 

analog  only 

4  Jun  1984 

19:16:48.9 

38.02 

87.03 

5.0 

2.3 

analog  only 

30  Jul  1984 

12:57:25.7 

38.16 

87.37 

10.0 

1.7 

analog  only 

1  Sep  1984 

18:27:44.2 

38.04 

87,29 

10.0 

1.8 

analog  only 

1  May  1985 

01:16:27.9 

37.99 

87.61 

1.4 

2.9 

digital  • 

10  Jan  1986 

19:54:52.0 

38.16 

87.58 

10.0 

2.5 

analog  only 

13  Feb  1987 

16:44:51.2 

38.18 

87.47 

5.0 

2.5 

digital 

16  Feb  1987 

16:27:06.2 

38.18 

87.22 

5.0 

2.5 

analog  only 

16  Jun  1992 

19:28:49.0 

38.06 

87.46 

2.0 

2.3 

digital  * 

16  Jun  1992 

23:11:13.6 

38.11 

87.39 

0.1 

2.65 

strip  mine  blast  • 

24  Jul  1992 

15:23  15.4 

38.11 

87.31 

0.1 

2.50 

strip  mine  blast  * 

*  relocated  with  CMVSN,  BLO,  and  KY  stations;  however  depths  are  from  archive  files 


Each  of  the  four  events  selected  for  further  processing  where  relocated 
as  part  of  this  study.  When  only  the  CMVSN  stations  were  used,  the  earth¬ 
quakes  relocated  with  depths  between  10-12  km  and  the  strip-mine  blasts  at 
0.0  km  (at  the  surface).  After  adding  BLO  and  the  Kentucky  stations,  it  was 
interesting  to  note  that  the  earthquakes  relocated  at  relatively  shallow 
depths  (=  2.0  km)  while  the  strip-mine  blasts  located  at  much  deeper  depths 
(6-12  km)  (Table  2).  The  velocity  model  for  this  area  will  need  to  be  improved 
to  accurately  determine  depths  to  the  earthquakes  without  portable  stations 
in  the  field. 

Waveforms  for  the  four  selected  events  are  shown  in  Figure  3  for  sta¬ 
tions  WDIN  and  BPIL.  Each  trace  is  preceded  by  1.0  second  of  background 
noise  prior  to  the  P-wave  arrival.  Strip-mine  blasts  occurring  in  the  Illinois 
Basin  are  distinctive  and  easily  recognizable  from  earthquakes  (Figure  2  and 
Figure  3).  All  have  a  low-frequency  Rg  phase  (a  fundamental  mode  surface 
wave  with  a  period  of  0.25  to  3.0  seconds),  usually  with  a  larger  amplitude 
Airy  phase.  Higher  mode  Rg  phases  are  occasionally  discernible  in  the  time 
series.  The  shape  of  the  Rg  phase  coda  varies  with  distance  and  source  loca¬ 
tion  (due  to  differences  in  sediment  thickness)  within  the  Illinois  Basin 
(Hutchenson  et  ai,  1990;  Hutchenson  and  Herrmann,  1991).  A  distinctive  Rg 
phase  and  coda  are  characteristic  of  shallow-source  events  (B&th,  1975; 
Shapiro,  1988;  Kafka,  1990).  Strip-mine  blasts  in  open-pit  mines  are  shallow 
source  events,  although  not  necessarily  point  sources. 
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Fig.  3.  A  sequence  of  four  time  histories  recorded  at  BFIL  and  WDIN.  For  each  station,  the 
upper  two  traces  correspond  to  two  different  blast,  and  the  lower  two  traces  to  two  different 
earthquake  recordings.  The  June  16,  1992  earthquake  and  blast  time  histories  are  included. 
Annotation  shows  distance  (DIST),  azimuth  (AZ),  maximum  ground  motion  (in  cm/sec)  and 
sample  rate  (DT).  The  low  frequency  Rg  coda  is  very  pronounced  in  the  blasts  at  WDIN.  The 
Rg  coda  is  not  as  pronounced  at  BPIL  for  the  two  blasts,  but  the  S  phase  arrivals  on  the  two 
earthquakes  are  more  distinct. 
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Earthquakes  generally  occur  at  greater  depths  than  blasts.  As  such, 
earthquakes  do  not  usually  have  the  characteristic  low-frequency  Rg  coda 
associated  with  shallow-source  events.  Reported  shallow- source  earthquakes 
are  discussed  later  in  this  section.  The  frequency  content  of  earthquakes  is 
higher,  with  P-,  S-,  and  sometimes  Lg,  as  the  only  identifiable  phases. 
Finally,  the  earthquake  codas  decay  much  more  rapidly. 

The  signals  of  the  earthquake  and  strip-mine  blast  of  16  June  1992 
were  first  examined  using  multiple  filter  analysis  to  permit  something  like  a 
sonogram  display  (Dziewonski  et  al.,  1969).  Herrmann  (1973,  1987)  modified 
the  technique  to  estimate  spectral  amplitudes  of  various  modes.  This  moving 
window  analysis  is  normally  used  for  obtaining  surface-wave  dispersion  but 
also  indicates  apparent  velocities  of  body  waves.  The  displacement  in  the 
time  domain  caused  by  a  dispersed,  propagating,  multimode  wave  at  a  dis¬ 
tance  r  can  be  represented  by 

oo 

fit, r)  =  2~  J  F(a>.  r)  e,a/t  d co  (1) 


or 


°°  M 

fit,  r)  =  -L  J  r)  e'*"1 "  kjf)  dco  (2) 

~oo 


where  t  is  time,  co  is  angular  frequency,  kj  is  the  wave  number  of  the  j’th 
mode,  Aj  is  the  complex  amplitude  of  the  j’th  mode,  and  M  is  the  highest 
mode  in  the  signal  (a  total  of  M  +  1  modes).  The  signal  is  then  filtered  with  a 
symmetric  Gaussian  band-pass  signal  H(<y)  of  the  form 

[exp(  -  a  (a)  -  (oQ)2  /  co*)  \co-co0^coc\ 

H(<y)  =  in  ,  ,  (3) 

0  I  co  -  co0  >  <yc  I 


where  a  is  a  filter  parameter,  coa  is  the  filter  center  frequency,  and  co  -  co0  ±  coc 
is  the  cutoff  frequency  for  the  filter.  The  resultant  time  filtered  signal  repre¬ 
senting  a  multi-mode  signal  is  of  the  form 


(on 


sl/2  M 

-  SAjl^o.r)- 
a )  j=o 


(4) 


exp[i(o>0t  -  kojr)}  exp 


where  U0J  is  the  group  velocity  of  the  j’th  mode  at  frequency  co  -  co0. 


-64- 


The  envelope  of  the  signal  is  formed  using  the  Hilbert  transform.  The 
group  velocity  of  a  particular  phase  is  determined  from  the  time  of  the  peak 
envelope  amplitude.  To  define  or  separate  various  phases  or  modes,  a  mini¬ 
mum  separation  of  2td  is  required,  where 


T0  is  the  period  corresponding  to  frequency  co0.  For  separations  less  than  2td, 
the  observed  envelope  maxima  will  not  be  the  true  mode  maxima,  but  a  com¬ 
plex  sum  of  the  modal  amplitudes  (Herrmann,  1973). 

The  envelope  maxima  for  both  the  earthquake  and  strip-mine  blast 
recorded  at  WDIN  on  16  June,  1992,  are  shown  in  Figures  4a  and  4b.  In  this 
display,  the  contours  of  constant  amplitude  are  shown.  The  small  square,  cir¬ 
cle  and  triangle  and  plus  symbols  indicate  the  largest  through  the  fourth 
largest  spectral  peaks,  respectively,  at  a  given  filter  period.  The  apparent  P- 
phase  velocity  of  the  earthquake  is  shown  with  a  velocity  between  6.0  to  6.3 
km/sec,  while  that  of  the  blast  is  slightly  slower,  a  high  frequency  arrival 
near  6.3  km/sec,  similar  to  the  earthquake,  but  with  strong  secondary 
arrivals  between  5.6  to  6.0  km/sec.  At  this  short  distance,  this  may  reflect  a 
slightly  different  ray  path  from  the  source  to  station  because  of  the  shallow¬ 
ness  of  the  blast  and  because  of  the  very  thick  sedimentary  rock  section  in 
the  region  (  >  3  km). 

The  Lg/S-phase  arrival  times  for  each  event  are  similar,  both  traveling 
near  3.5  to  3.75  km/sec.  However,  note  the  large  difference  in  frequency  con¬ 
tent  in  the  signal.  The  Lg-phase  arrival  of  the  earthquake  is  not  obvious 
below  15  Hz  while  the  strip-mine  blast  shows  several  maxima  between  3  to 
15  Hz  with  little  signal  content  above  15  Hz. 

Two  body  wave  arrivals  are  shown  on  both  the  earthquake  and  strip- 
mine  blast  with  velocities  near  4.6  and  5.0  km/sec.  These  arrivals  are  also 
observable  in  the  time  series  waveform.  Both  arrivals  exhibit  similar  veloci¬ 
ties,  regardless  of  source  location.  The  arrivals  may  be  reflections  or  refrac¬ 
tions  from  several  of  the  known  groups  in  the  vicinity,  for  example,  the  Chat¬ 
tanooga  shale  or  Eau  Claire  formation  (Sexton  et  al.,  1986). 

The  fundamental  mode  Rg-phase  and  at  least  one,  strong  higher  mode 
are  visible  for  the  blast.  The  fundamental  mode  is  visible  (extreme  lower  left- 
hand  corner  of  the  contour  plot)  with  a  group  velocity  between  1.0  to  2.0 
km/sec  and  frequency  between  0.5  and  2.0  Hz.  The  higher  mode  is  visible 
(the  next  maxima  towards  the  upper  right  of  the  figure)  with  a  group  velocity 
between  1.5  to  2.9  km/sec  and  frequency  of  1.0  to  3.2  Hz. 

Other  contours  represent  incoherent  energy  due  to  either  noise  or  scat¬ 
tering.  As  a  particular  feature  of  the  contouring,  the  incoherent  energy  may 
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Fig.  4a.  Envelope  maxima  as  a  function  of  frequency  for  the  earthquake  of  16  June  1992  at 
WDIN.  The  40.96  seconds  of  recorded  signal  are  shown  on  the  far  right  in  linear  time,  lb 
its  left  is  the  same  signal  plotted  as  a  function  of  group  velocity  to  assist  in  relating  the 
observed  signal  arrivals  to  specific  group  velocities. 

be  at  a  higher  level  than  any  coherent  energy.  For  example,  note  the  low  fre¬ 
quency  (0.70  to  3.0  Hz)  band  on  the  earthquake  (Figure  4a).  It  is  not  coher¬ 
ent  but  is  low  frequency  information  with  an  apparent  group  velocity  greater 
than  2.0  km/sec.  Likewise,  the  points  along  the  top  margin  at  all  frequencies 
represent  noise  energy  with  apparent  group  velocities  greater  than  6.5 
km/sec.  The  apparent  absence  of  information  in  the  3  to  12  Hz  band  is  a 
direct  result  of  the  incoherent  energy  mapped  along  the  top  of  the  figure.  In 
this  same  band  for  the  explosion  (Figure  4b),  the  coherent  information  is  at  a 


Fig.  4b.  Envelope  maxima  as  a  function  of  frequency  for  the  strip-mine  blast  of  16  June  1992 
at  WDIN.  A  well-developed  fundamental  and  at  least  one  higher  mode  are  visible  in  the 
lower  left  comer  of  the  group  velocity  contours,  between  1.0  to  2.9  km/sec  and  0.5  to  3.2  Hz. 


higher  level  than  the  incoherent  noise.  Information  observable  from  the  con¬ 
tours  is  also  directly  observable  in  individual  spectra. 

The  spectral  content  of  the  P,  Lg,  and  Rg  phases  were  examined  in  addi¬ 
tional  detail  at  both  WDIN  and  BPIL.  Each  spectra  contains  512  data  sam¬ 
ples,  obtained  from  a  velocity  window  using  the  limits  defined  by  the  group 
velocity  maxima  of  each  phase  (Dziewonski  et  al.,  1959;  Herrmann,  1973). 
The  P-phase  spectra  were  windowed  between  5.0  and  6.5  km/sec  at  both  sta¬ 
tions.  The  Lg  phase  spectra  were  windowed  between  2.5  to  4.0  km/sec  at 
WDIN  and  3.0  to  4.0  km/sec  at  BPIL.  Finally,  the  Rg  phase  spectra  were 
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windowed  between  0.5  to  2.5  km/sec,  the  window  nearly  overlapping  with 
those  of  the  Lg  phase  at  WDIN. 

The  windowed  time  series  were  detrended,  then  tapered  with  a  10% 
cosine  taper  applied  to  each  end  of  the  window  to  help  remove  high  frequen¬ 
cies  associated  with  edge  effects.  A  Fast  Fourier  Transform  technique  (Her¬ 
rmann,  1987)  was  used  to  calculate  each  of  the  displacement  spectra.  Finally, 
the  spectra  were  corrected  to  ground  motion  at  the  observed  distance  by 
removing  instrument  effects.  No  distance  correction  is  applied  since  source 
estimates  are  not  made.  Spectra  for  each  phase,  event,  and  station  are 
shown  on  the  same  scale  in  Figures  5  through  8. 


Fig.  5.  P-phase  spectra,  WDIN  for  the  traces  shown  in  Figure  3. 
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The  spectra  for  each  phase  and  station  confirm  the  frequency  content 
observed  in  the  seismograms  (Figures  2  and  3).  These  earthquakes  generally 
have  P-phase  comer  frequencies  near  10  Hz,  with  significant  energy  up  to 
20-30  Hz.  Strip-mine  blasts,  however,  have  much  lower  P-phase  comer  fre¬ 
quencies,  closer  to  4-6  Hz,  and  high  frequency  energy  is  usually  not  obvious 
above  10  Hz. 

A  comparison  of  the  recordings  of  the  same  event  at  the  different  sta¬ 
tions  shows  similarity  in  spectra,  with  the  exception  of  the  920724  blast 
because  of  the  low  S/N  ratio  at  BPIL.  In  addition,  at  a  given  station,  the  spec¬ 
tra  of  the  two  blasts  and  of  the  two  earthquakes  are  also  very  similar.  Thus 
the  site  effects  at  the  stations  are  either  not  affecting  the  spectral  shapes  or 
are  affecting  them  in  the  same  manner.  The  blast  spectra  at  BPIL  and  WDIN 
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Fig.  7.  Lg-phase  spectra,  WDIN  for  the  traces  shown  in  Figure  3. 


show  a  pronounced,  isolated  spectral  peak  btween  1.5  and  3.5  Hz.  Since  this 
is  not  as  apparent  in  the  noise  or  earthquake  spectra,  this  must  be  a  feature 
of  the  strip-mine  explosion.  In  fact,  the  modulation  of  the  WDIN  recordings  of 
the  920616  blast  between  1  and  8  Hz  is  characteristic  of  that  due  to  ripple  fir¬ 
ing  (Baumgardt  and  Ziegler,  1988). 

The  significant  conclusion,  then,  from  examining  the  P-wave  spectra  is 
that  the  earthquakes  are  richer  in  higher  frequencies  than  the  strip-mine 
blasts  with  higher  corner  frequencies.  This  is  opposite  of  the  results  one 
would  expect  from  theoretical  spectral  estimates  between  earthquakes  and 
nuclear  explosions  (Evemden  et  al.,  1986).  Studies  for  nuclear  explosions 
and  earthquakes  in  the  western  United  States  suggest  similar  findings  (Mur¬ 
phy  and  Bennett,  1982;  Bennett  and  Murphy,  1986,  Taylor  et  al,  1989;  Chael, 
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Fig.  8.  Lg-phase  spectra,  BPIL  for  the  traces  shown  in  Figure  3. 


1988, 1991). 

Results  for  the  Lg  phase  are  similar.  The  Lg  spectra  for  the  blasts  show 
significantly  higher  low  frequency  energy  (0.8  to  5  Hz)  while  the  earthquake 
spectra  are  richer  in  energy  at  higher  frequencies  (>  8  Hz).  This  is  similar  to 
the  trends  observed  in  the  P-wave  spectra.  A  spectral  peak  is  seen  in  some  of 
the  blast  spectra  in  the  1  -  3  Hz  range,  but  the  spectral  modulation  is  not 
readily  apparent.  Another  interesting  feature  is  that  the  Lg-  and  P-wave 
spectral  levels  are  about  the  same  for  the  same  station-event  tuple. 

The  Rg  and  Lg  spectral  results  are  similar  in  terms  of  frequency  con¬ 
tent.  However,  note  that  at  station  BPIL,  the  farthest  from  the  source,  the 
spectral  level  is  almost  indistinguishable  from  the  background  noise  for  the 
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Fig.  9.  Rg-phase  spectra,  WDIN  for  the  traces  shown  in  Figure  3. 


earthquake;  only  the  strip-mine  blast  spectra  shows  energy  in  the  0.4  to  2  Hz 
band. 

No  focal  mechanisms  have  been  obtained  for  earthquakes  in  this  area 
due  to  the  size  of  the  events  and  the  CIvTVSN  station  distribution.  It  is 
assumed  in  this  study  that  the  earthquakes  all  have  similar  mechanisms, 
with  fault  planes  trending  NNE  or  NNW.  This  assumption  is  based  on  the 
general  NNE  trend  of  mapped  faults  and  basement  structures  in  the  immedi¬ 
ate  area  (Pratt  et  al.,  1989;  Sexton  et  al.,  1986;  Nelson  and  Lumm,  1984),  or 
from  extending  the  La  Salle  Anticlinal  Belt  through  the  basement  in  this 
area  (Hamburger  and  Rupp,  1988).  Most  workers  agree  that  current  seismic¬ 
ity  is  the  result  of  reactivated  N-S  trending  basement  faults,  dissimilar  to  the 
E-W  structures  of  the  Cottage  Grove-Rough  Creek  Fault  Systems  (Nelson 
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Fig.  10.  Rg-phase  spectra,  BPIL  for  the  traces  shown  in  Figure  3. 


and  Lumm,  1984).  The  few  focal  mechanisms  for  the  nearest  earthquakes  in 
Illinois  (Herrmann,  1979b;  Taylor  et  al. ,  1989)  may  be  associated  with  the  La 
Salle  Anticlinal  Belt  (Hamburger  and  Rupp,  1988).  Thus,  the  relative  event 
levels  at  a  given  station  should  reflect  effects  of  the  source-time  function,  not 
radiation  patterns.  If  true,  then  the  earthquake  spectra  can  be  directly  com¬ 
pared. 

This  may  not  be  the  case  for  the  strip-mine  blasts.  Many  of  the  strip 
mine  shots  are  patterned  averaging  three  rows,  two  decks,  and  approxi¬ 
mately  a  hundred  holes  per  shot.  Azimuthal  effects  may  be  important,  but 
they  cannot  be  resolved  with  this  data  set. 


-73- 


A  remaining  question  of  importance  is  whether  shallow  earthquakes  in 
the  region  exhibit  similar  Rg  codas.  One  hundred  one  (101)  events  were 
found  in  the  CMVSN  event  catalog  with  depths  less  than  4.0  km  between  37° 
to  40°  N  latitude  and  86.5°  to  92.5°  W  longitude.  Many  of  these  are  within 
the  Illinois  Basin  (Figure  11).  Well-located  earthquakes  in  the  basin  gener¬ 
ally  locate  at  average  depths  of  approximately  10.0  to  11.0  km.  The  catalog 
shows  a  strong  correlation  of  depths  at  2.5,  5.0,  and  10.0  km,  the  default  loca¬ 
tion  depths.  No  shallow-depth  earthquakes  have  been  reported  by  the 
CMVSN  in  Indiana;  however,  it  is  recognized  that  depth  is  poorly  constrained 
outside  the  coverage  area  of  the  CMVSN  stations.  All  events  used  in  this 
study  were  relocated,  with  particular  emphasis  on  depth. 


Fig.  11.  Shallow  earthquakes  reported  in  the  the  Illinois  Basin  and  eastern  Missouri  in  rela¬ 
tion  the  CMVSN  stations. 

Waveforms  for  several  reported  "shallow-depth"  events  (Table  3)  are 
shown  in  Figures  12  through  14  at  various  CMVSN  stations.  Although  nomi¬ 
nally  shallow  (<  4.0  km),  none  have  the  low-frequency  Rg  coda  typical  of  shal¬ 
low  events.  The  implication  is  that  either  the  sources  were  deeper  than  the 
layered  sedimentary  rock  waveguide,  or,  more  likely,  the  published  bulletin 
depths  are  in  error  due  to  poor  depth  control. 

Shallow  earthquakes  do  occur  within  or  on  the  boundaries  of  the  Illi¬ 
nois  Basin.  An  earthquake  that  occurred  on  14  August  1965  in  southern  Illi¬ 
nois  produced  a  short  period  dispersed  surface  wave  on  a  World  Wide 
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Fig.  12.  Time  series  of  the  earthquake  of  29  June  1984  in  southeastern  Illinois  at  various 
CMVSN  stations.  The  traces  are  annotated  with  station  name,  peak  velocity  in  cm/sec,  epi- 
central  distance,  and  azimuth  (°’s  relative  to  north).  Note  that  the  waveforms  at  stations 
FVMZ  and  SLMZ  have  a  noticeably  different  character.  Both  have  a  narrower  instrument 
response  with  a  lower  passband.  As  a  result,  the  high  frequency  information  seen  in  the 
other  CMVSN  stations  is  not  observed.  However,  short  period  Rg  waves  are  not  seen  either. 

Standard  Seismic  Network  (WWSSN)  long-period  (LP)  instrument  (FLO  in 
St.  Louis)  for  a  path  length  of  200  km  (Herrmann,  1974;  Herrmann  and  Nut- 
tli,  1975).  This  observation  was  possible  only  because  the  LP  instrument  fil¬ 
tered  out  the  high-frequency  Lg  signal.  This  mb  3.8  event  was  assigned  a 
depth  of  1.0  km  based  on  surface-wave  modeling  (Herrmann,  1974).  No 
broader  bandwidth  recordings  of  this  event  exist  since  the  CMVSN  was  not 
installed  until  1974  (Stauder  et  al.,  1991). 
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Fig.  13.  Time  series  of  the  earthquake  of  24  October  1990  in  centra!  Illinois  at  various 
CMVSN  stations. 

DISCUSSION  AND  CONCLUSIONS 

The  occurrence  of  earthquakes  in  the  same  region  as  active  strip  mines 
allows  a  direct  comparison  of  the  appearance  and  spectral  content  of  each 
event  type.  Arguments  that  differences  in  signal  character  of  body  waves  are 
due  to  path  propagation  effects  are  weakened  when  both  types  of  events  are 
located  in  the  same  vicinity,  assuming  that  the  depth  effect  is  small,  and 
when  the  same  stations  record  both  types  of  events.  Events  similarly  located 
also  reduce  the  variability  found  in  phase  onsets,  phase  appearance,  and  par¬ 
ticular  phases  found  from  various  sites  (Gupta  and  Hartenberger,  1981). 
Results  in  this  study  are  similar  in  these  respects  to  similar  studies  using 


nuclear  explosion  and  earthquakes  (i.e.  Chad,  1988,  1991). 

On  analog  or  digital  recordings  (Figures  2  and  3),  the  appearance  of  the 
strip-mine  blast  and  earthquake  signals  allow  for  visual  discrimination 
between  the  two.  Signals  from  strip-mine  blasts  have  a  distinctive  Rg  phase, 
most  with  a  large-amplitude  Airy  phase,  indicative  of  shallow-source  events 
(B$th,  1975;  Shapiro,  1988;  Kafka,  1990).  Earthquakes,  even  presumably 
shallow  earthquakes  in  the  Illinois  Basin,  lack  the  well-developed  Rg  phase 
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and  have  rapidly  decaying  codas.  The  lack  of  a  well-developed  Rg  phase  in 
reported  shallow  earthquakes  quite  probably  indicates  poor  depth  control  in 
the  earthquake  catalog;  these  events,  reported  at  a  shallow  depth,  are  thus 
inferred  to  occur  much  deeper.  At  least  one  demonstratably  shallow  earth¬ 
quake  prior  to  the  installation  of  the  CMVSN  has  produced  surface  waves 
(Herrmann,  1974;  Herrmann  and  Nuttli,  1975),  indicating  most  earthquakes 
are  deeper.  Finally,  the  onset  of  the  P,  S,  and  sometimes  Lg,  phases  are  more 
impulsive  for  earthquakes  than  for  the  strip-mine  blasts. 

Aside  from  visual  discrimination,  the  spectral  content  also  shows  differ¬ 
ences.  The  corner  frequencies  of  the  P-wave  spectra  are  lower  for  the  blast 
than  for  the  earthquakes.  In  addition,  the  energy  content  in  the  high- 
frequency  passband  (10-30  Hz)  is  higher  for  earthquakes  than  for  strip-mine 
blasts.  This  result  is  similar  to  prior  studies  of  regional  discriminants  in  the 
Western  U.S.  and  Central  Asia  (Balapan)  designed  to  distinguish  between 
nuclear  explosions  and  earthquakes  (Bennett  et  al.,  1992).  These  studies 
have  shown  both  Pg  and  Lg  spectra  richer  in  higher  frequency  energy,  espe¬ 
cially  between  3  to  6  Hz  (Murphy  and  Bennett,  1982;  Bennett  and  Murphy, 
1986;  Bennett  et  al,  1992),  in  comparisons  of  spectral  ratio  in  the  1  to  2  Hz 
and  6  to  8  Hz  passbands  (Taylor  al.,  1989),  and  at  frequencies  between  10 
to  30  Hz  (Chael,  1988,  1991). 

Results  for  Lop  Nor,  the  Chinese  test  site,  indicate  the  opposite  (Bennett 
et  al.,  1992).  Earthquakes  exhibit  greater  low  frequency  energy  while  the 
nuclear  explosions  were  enhanced  at  the  higher  frequencies.  In  addition,  the 
raw  waveforms  showed  a  more  prominent  Rg  phase  (larger  amplitude,  but 
not  larger  than  S  or  Lg)  than  the  explosions  at  one  of  the  stations. 

Theoretical  source  models  suggest  that  explosions  should  be  richer  in 
high  frequency  energy  primarily  due  to  the  differences  in  corner  frequency 
and  spectral  rolloff  rates  (Evemden  et  al.,  1986;  Walter  and  Priestly,  1991). 
In  this  case,  the  earthquake  has  lower  high  frequency  content  than  the  explo¬ 
sion  at  fixed  low  frequency  levels  because  of  the  larger  source  dimension 
(lower  stress  drop)  of  the  earthquake  relative  to  the  point  source  nuclear 
explosion.  Our  study  finds  the  opposite  effect  for  the  same  reason.  For  a 
given  high  frequency  level,  used  to  assign  the  event  magnitude,  the  earth¬ 
quake  appears  as  a  point  source  while  the  strip-mine  explosion  appears  as  an 
extended  source.  This  is  not  surprising,  given  the  large  spatial  dimensions  of 
the  strip-mine  shots  as  well  as  the  time  delays  inherent  to  ripple-firing. 

This  observation  may  serve  to  distinguish  point  sources  from  dis¬ 
tributed  sources  for  small  events.  The  implication  in  a  proliferation  environ¬ 
ment  is  that  one  must  look  at  attributes  other  than  P-  or  Lg-wave  spectral 
content  to  define  the  nuclear  explosion  (path  attenuation,  earthquake  radia¬ 
tion  patterns,  etc.).  The  appearance  of  the  Rg  wave  and  signal  complexity 
may  provide  an  answer 
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Finally,  the  occurrence  of  earthquakes  in  this  region,  and  elsewhere 
within  the  Illinois  Basin,  are  apparently  the  result  of  movement  along  exist¬ 
ing  Paleozoic  basement  structures  (Hamburger  and  Rupp,  1988;  Pratt  et  al., 
1989;  Taylor,  1991).  While  a  concern  of  the  public,  it  cannot  be  reasonably 
argued  that  earthquakes  in  this  region  are  a  result  of  any  current  blasting  or 
mining  activity.  Unlike  other  areas  where  a  spatial  association  with  mining 
exists  (ie  McGarr,  1984;  Gibowicz,  1990),  no  temporal  association  between 
earthquakes  and  strip-mine  blasts  have  been  found  in  this  area.  The  spatial 
association  is  due  to  the  location  of  Paleozoic  basement  structures  beneath 
easily  obtainable  surface  coal  resources. 

While  seismicity  has  been  associated  with  active  mining  operations 
(McGarr,  1984;  Gibowicz,  1990),  it  is  rare  with  surface  mining  (Gibowicz, 
1990).  Events  have  occurred  as  a  result  of  coal  mining  operations  in  Poland. 
However,  the  induced  events  in  that  area  appear  to  be  triggered  by  two 
actions,  a  decrease  in  the  vertical  stress,  caused  by  the  removal  of  the  over¬ 
burden  (nearly  300  meters)  and  an  increase  in  the  effective  stress,  due  to 
decreasing  pore  pressure  as  a  result  of  groundwater  withdrawal  (Gibowicz, 
1990).  Obtainable  coal  resources  (strip-mines,  not  underground  mines) 
occurs  at  shallower  depths  in  Illinois  and  Indiana,  about  a  third  less  overbur¬ 
den  is  removed. 
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DISCOURSE  ON  THE  USE  OF  SHORT  PERIOD  SURFACE  WAVES 
AS  A  DEPTH  DISCRIMINANT 


R.B.  Herrmann 
ABSTRACT 


Consideration  of  the  usefulness  of  short-period  surface  wave 
recordings  in  a  regional  discrimination  environment  is  pre¬ 
sented.  These  signals  are  not  evidence  of  an  explosion,  but 
rather  of  shallow  source  depth.  If  independent  estimates  of  the 
event  seismic  moment  are  available,  the  presence  or  absence  of  a 
short-period  recording  can  be  used  as  a  depth  constraint  on  the 
event,  which  of  itself  may  be  a  strong  discriminant. 

INTRODUCTION 

It  is  accepted  that  the  observation  of  a  short  period  surface  wave  is  am 
indication  of  a  shallow  source  depth.  This  general  statement  is  based  on 
numerous  observations  of  mining  activity  which  produce  the  characteristic 
short  period  surface  waves  seen  on  seismograms.  The  theoretical  basis  lies  in 
the  fact  that  the  surface  wave  excitation  can  be  related  to  the  product  of  an 
eigenfunction  sampled  at  source  depth  and  another  at  receiver  depth  divided 
by  a  depth  integral  of  an  eigenfunction  squared.  Since  short  period  surface 
wave  eigenfunctions  do  not  penetrate  very  deep,  the  surface  wave  recorded  at 
the  surface  and  generated  by  a  shallow  source  will  be  significantly  larger 
than  one  generated  by  a  deeper  source  (Tsai  and  Aki,  1969;  Herrmann,  1974). 

This  observation  is  typically  used  by  regional  seismologists  to  classify 
events  as  not  being  earthquakes.  The  strong  assumption  is  that  shallow 
earthquakes  either  do  not  occur,  or  have  sufficient  signal  complexity  to  over¬ 
whelm  the  coherent  short  period  surface  wave.  This  is  not  always  true. 

Figure  1  presents  WWSSN  long-period  recordings  of  a  small,  mbLg  =  3.8, 
earthquake  that  occurred  near  Cairo,  Illinois  on  August  14,  1965  (Herrmann 
and  Nuttli,  1975).  This  earthquake  was  only  felt  over  a  short  distance  of  15 
km,  but  did  cause  intensity  VII  damage  (Nuttli  and  Zollweg,  1974)  and  led  to 
Nuttli’s  (1973)  development  of  the  magnitude  scale.  This  was  an  earth¬ 
quake  since  no  large  coal  strip  mines  were  nearby  which  would  have  shot  off 
greater  than  100  tons  of  explosive.  This  earthquake  did  generate  observable 
surface  waves  on  long  period  WWSSN  seismographs  1400  km  away  (Her¬ 
rmann,  1974)  and  analysis  yielded  strike-slip  focal  mechanism,  seismic 
moment  of  2. 9 1021  dyne-cm  and  a  depth  of  1.5  km  (Herrmann,  1979). 
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OBSERVED 


Fig.  1.  Comparison  of  observed  and  predicted  waveforms  recorded  on  the  OXF  and  FLO 
long-period  instruments  for  the  August  14,  1965  earthquake.  The  stations  are  300  and  200 
km,  respectively,  from  the  source.  The  15-100  WWSSN  LP  seismographs  had  gains  of  3. OK. 

We  conclude  that  shallow  earthquakes  can  also  generate  short  period 
surface  waves  that  can  be  observed  at  regional  distances.  Hence  one  must  be 
cautious  about  their  use  as  an  event  type  discriminant. 

Another  statement  made  concerning  regional  recordings  is  that  there  is 
little  value  in  studying  the  short-period  surface  wave  since  it  does  not  propa¬ 
gate  to  large  distances.  The  previous  example  implies  that  under  proper  con¬ 
ditions  of  correct  recording  instrument  and  band  and  with  a  relatively  uni¬ 
form  waveguide,  as  in  the  central  United  States,  these  waves  can  be 
observed. 

If  this  is  the  case,  then  can  they  be  useful  for  the  discrimination  prob¬ 
lem? 

DEPTH  CONSTRAINTS 

The  key  point  is  that  the  presence  short-period  surface  wave  is  a  source 
depth  indicator  and  not  a  source  classifier.  In  a  discrimination  context,  if  it 
can  be  shown  that  an  event  is  greater  than  3  km,  perhaps,  then  the  event  is 
most  likely  not  an  explosion.  The  event  can  then  be  confidently  identified  as 
an  earthquake  and  eliminated  from  further  consideration.  In  the  today’s 
environment  of  automatic  event  location  and  classification,  it  would  be  inter¬ 
esting  to  consider  how  surface  wave  information  can  be  incorporated. 

The  first  step  is  to  generate  synthetic  surface-wave  seismograms  for 
point  earthquake  and  explosion  sources  as  a  function  of  source  depth.  The 
ground  motions  generated  can  then  be  passed  through  different  instruments 
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Synthetic  seismograms  were  generated  at  distances  of  100  and  1000 
km  for  a  point  explosion  and  for  a  strike-slip  earthquake  source  observed  at 
an  azimuth  of  22°  from  strike.  The  source  time  function  had  a  duration  of  0.4 
seconds,  which  would  not  affect  the  long  period  recordings.  It  also  should  not 
affect  the  1  Hz  surface  wave  on  the  short  period  instrument.  The  duration 
may  also  be  typical  of  an  event  with  seismic  moment  of  10  a  dyne-cm. 

To  present  the  results  of  the  computations  in  a  form  useful  for  discrimi¬ 
nation  studies,  a  plot  is  made  of  the  seismic  moment  required  to  produce  a 
threshold  amplitude  on  the  seismogram  for  the  event  depth.  This  threshold 
is  taken  to  be  0.1  cm,  which  is  the  limit  of  what  can  be  seen  and  measured  on 
photographic  seismograms.  Figure  2  presents  such  information.  Figure  2a  is 
for  the  long-period  instrument  and  Figure  2b  is  for  the  short  period  instru¬ 
ment.  The  striking  feature  is  that  the  explosion  and  strike-slip  earthquake 
give  identical  observed  amplitudes  at  shallow  depth.  Beneath  the  Paleozic 
layer,  differences  are  seen,  with  the  explosion  generating  smaller  surface 
waves. 

These  figures  can  be  used  very  simply  if  the  event’s  seismic  moment  is 
known.  If  the  observed  amplitude  is  greater  than  the  threshold,  then  its 
depth  must  be  to  the  left  of  the  curve;  if  the  signal  is  not  observed,  then  the 
depth  is  in  the  region  to  the  right  of  the  curve.  This  is  a  simple  binary  deci¬ 
sion.  Changing  the  threshold  will  only  change  the  set  of  master  cuves.  To 
account  for  unknowns  in  the  earth  model,  a  suite  of  curves  could  be 


Fig.  2.  (a)  Plot  of  moment  required  to  exceed  a  recorded  surface-wave  amplitude  of  0.1  cm  on 
a  1.5K  15-100  vertical  component  WWSSN  LP  seismograph  at  distances  of  100  and  1000 
km,  and  for  explosion,  EX,  and  vertical  strike-slip,  SS,  sources.  The  receiver  azimuth  from 
the  strike-slip  source  is  22  , ‘3°  from  strike. 

(b)  Plot  of  moment  required  to  exceed  a  recorded  surface- wave  amplitude  of  0.1  cm  on  a  50K 
vertical  component  WWSSN  SP  seismograph  at  a  distance  of  100  km,  and  for  explosion,  EX, 
and  vertical  strike-slip,  SS,  sources.  The  receiver  azimuth  from  the  strike-slip  source  is  22.5° 
from  strike. 

generated.  This  would  lead  to  a  distribution  in  possible  depth  for  a  given  seis¬ 
mic  moment.  The  absence  of  a  signal  at  the  threshold,  perhaps  due  to  its 
being  obscured  by  noise,  could  then  be  assigned  a  conditional  probability  that 
its  depth  is  greater  than  a  certain  value.  The  same  distribution  could  be  used 
with  an  actual  observed  amplitude  to  assign  a  depth  through  another  set  of 
curves.  When  combined  with  some  a  priori  statement  that  an  event  with 
depth  greater  than  3  km,  for  example,  is  not  an  explosion,  then  the  result 
would  be  a  likelihood  that  the  event  is  not  an  explosion. 

DISCUSSION 

This  simple  use  of  master  curves  shows  how  short  period  surface  wave 
observations  could  be  used  to  quantify  depth  estimates.  The  requirements  are 
as  follow: 

a)  the  regional  surface-wave  propagation  characteristics  are  known,  e.g., 

the  crustal  velocity  and  Q  model  is  known. 
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b)  there  is  an  independent  estimate  of  the  seismic  source  size,  e  g.,  a 

moment-magnitude  relationship. 

The  construction  of  master  curves  is  simple  with  existing  code.  The  incorpo¬ 
ration  of  this  technique  into  a  rules  based  depth  classification  process  is  not 
too  difficult. 

Figure  2a  indicates  that  it  will  be  difficult  to  distinguish  between  shal¬ 
low  explosions  and  earthquakes  on  the  basis  of  long  period  recordings.  This 
permits  an  indirect  thought  experiment  to  tie  an  explosion  yield  to  an  mLg 
value  for  a  small  event.  Herrmann  et  al  (1992)  analyzed  a  large  number  of  1 
ton  point  chemical  explosions  and  found  that  the  isotropic  moment  was 
roughly  21018  dyne-cm.  Denny  and  Johnson  (1991)  proposed  a  one-to-one 
relation  between  seismic  moment  and  yield.  Thus  the  seismic  moment  for 
the  August  14,  1965  event  would  be  indicative  of  a  1  kT  event.  It  is  interest¬ 
ing,  coincidental  and  certainly  not  well  understood,  that  the  assigned  mbLg  = 
3.8  number  is  roughly  within  the  proper  range  for  a  1  kT  nuclear  explosion 
Thus  one  should  be  able  to  observe  not  only  Lg  but  also  surface  waves  from 
such  a  nuclear  event  in  at  regional  distances. 

There  is  another  aspect  to  1  kT  explosions.  They  should  be  well 
observed  in  the  4  -  40  second  period  band  if  the  upper  crust  is  fairly  uniform. 
Thus  an  examination  of  radiation  patterns  of  the  Love  and  Rayleigh  waves 
would  be  able  to  distinguish  the  explosion  from  the  earthquake.  Thus  may  be 
the  threshold  for  doing  this  because  of  reduced  S/N  for  smaller  events  due  to 
the  presence  of  microseisms. 


REFERENCES 

Denny,  M.  D.,  and  L.  R.  Johnson  (1991).  The  explosion  seismic  source  func¬ 
tion:  models  and  scaling  laws  reviewed,  in  Explosion  Source  Phenomenol¬ 
ogy,  S.  R.  Taylor,  H.  J.  Patton  and  P.  G.  Richards  (eds),  American  Geophys¬ 
ical  Union,  Washington,  268  pp, 

Herrmann,  R.  B.  (1974).  Surface  wave  generation  by  central  united  States 
earthquakes,  Ph.  D.  Dissert.,  Saint  Louis  University,  263  pp. 

Herrmann,  R.  B.  (1979).  Surface  wave  focal  mechanisms  for  eastern  North 
American  earthquakes  with  tectonic  implications,  J.  Geophys.  Res.  84, 
3543-3552. 

Herrmann,  R.  B.,  and  0.  W.  Nuttli  (1975).  Ground  motion  modeling  in  a  con¬ 
tinental  interior,  I.  Theory  and  observations.  International  Journal  of 
Earthquake  Engineering  and  Structural  Dynamics  4,  49-58. 


-89- 


Herrmann,  R.  B.,  G.  Al-Eqabi,  and  K.  Hutchensen  (1992).  Isotropic  moment- 
yield  relations  for  small  chemical  explosions,  in  Quantification  of  mt ,  for 
small  explosions.  Scientific  Report  No.  1,  PL-TR-92-2109,  Phillips  Labora¬ 
tory,  Air  Force  Systems  Command,  Hanscom  Air  Force  Base,  Mas¬ 
sachusetts,  47  pp.  ADA253915 

Nuttli,  O.  W.  (1973).  Seismic  wave  attenuation  and  magnitude  relations  for 
eastern  North  America,  J.  Geophys.  Res.  78,  876-885. 

Nuttli,  O.  W.,  and  J.  E.  Zollweg  (1974).  The  relation  between  felt  area  and 
magnitude  for  central  United  States  earthquakes,  Bull.  Seism.  Soc.  Am. 
64,  73-85. 


Tsai,  Y.  B.,  and  K.  Aki  (1970).  Precise  focal  depth  determination  from  ampli¬ 
tude  spectra  of  surface  waves,  J.  Geophys.  Res.  75,  5729-5743. 


-90- 


Prof.  Thomas  Ahrens 

Seismoiogical  Lab,  252-21 

Division  of  Geological  &  Planetary  Sciences 

California  Institute  of  Technology 

Pasadena,  CA  91125 

Prof.  Kciiti  Aki 

Center  for  Earth  Sciences 

University  of  Southern  California 

University  Park 

Los  Angeles,  CA  90089-0741 

Prof.  Shelton  Alexander 
Geosciences  Department 
403  IV ike  Building 
The  Pennsylvania  State  University 
University  Park,  PA  19802 

Dr.  Ralph  Alewine,  111 
DARPA/NMRO 
3701  North  Fairfax  Drive 
Arlington,  VA  22203-1714 


Prof.  Charles  B.  Archambean 
ORES 

University  of  Colorado 
Boulder.  CO  80309 


Dr.  Thomas  C.  Bache,  Jr. 
Science  Applications  lnt'1  Corp. 
10260  Campus  Point  Drive 
San  Diego,  CA  92121  (2  copies) 


Prof.  Muawia  Barazangi 

Institute  for  the  Study  of  the  Continent 

Cornell  University 

Ithaca,  NY  14853 


Dr.  Jeff  Barker 

Department  of  Geological  Sciences 
Suite  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 

Dr.  Douglas  R.  Baumgardt 
EN SCO.  Inc 
54(H)  Port  Royal  Rond 
Springfield,  VA  22151-2388 


Dr.  Susan  Beck 
Department  of  Geosciences 
Building  #77 
University  of  Arizona 
Tuscon.AZ  85721 


Dr.  TJ.  Bennett 

S  (  CUM) 

A  Div  ision  of  Max  well  Laboratories 
1 1 800  Sunrise  Valiev  Drive  Suite  i?l 
Reston.  VA  22001  ' 

Dr.  Robert  Blandford 
AFTAC/TL  Center  for  Seismic  Suidie 
1 3(H)  North  1 7th  Street 
Suite  1450 

Arlington.  VA  22209-2308 

Dr.  Stephen  Bratt 
Center  for  Seismic  Studies 
1300  North  17th  Street 
Suite  1450 

Arlington.  VA  22209-2308 

Dr.  Lawrence  Burdick 
IGPP.  A -025 

Scripps  Institute  of  Oceanography 
Universitv  of  California,  San  Dieco 
La  Jolla.  CA  92093 

Dr.  Robert  Burridge 
Schlumbcrgcr  Doll  Research  Center 
Old  Quart  v  Road 
Ridgefield.  CT  06877 


Dr.  Jerry  Carter 
Center  for  Seismic  Studies 
1300  North  17th  Street 
Suite  1450 

Arlington.  VA  22209-2308 

Dr.  Eric  Chael 
Division  924! 

Sandia  Laboratory 
Albuquerque,  N\1  87185 


Dr.  Martin  Chapman 
Department  of  Geological  Sciences 
Virginia  Polytechnical  Institute 
21044  Derring  Hall 
Blacksburg,  VA  24061 

Prof.  Vernon  F.  Cormier 
Department  ofGeologv  &  Geophysics 
U-45,  Room  207 
University  of  Connecticut 
Storrs,  CT  06268 

Prof.  Steven  Day 

Department  of  Geological  Science-; 

San  Die  so  State  Universitv 
San  Diego,  CA  92182 


Marvin  Denny 
U.S.  Department  of  Energy 
Office  of  Amis  Control 
Washington,  DC  20585 


Dr.  Zoltan  Der 
ENSCO.  Inc. 

5400  Port  Roval  Road 
Springfield,  V A  22151-2388 


Prof.  Adam  Dziewonski 
Hoffman  Laboratory,  Harvard  University 
Dept,  of  Earth  Atmos.  &  Planetary  Sciences 
20  Oxford  Street 
Cambridge.  MA  02138 

Prof.  Jolm  F-bel 

Department  of  Geology  &  Geophysics 
Boston  College 
Chestmt  Hilf  MA  02167 


Eric  Fielding 
SNEE  Hall 
IN'STOC 

Cornell  University 
Ithaca,  NY  14853 

Dr.  Mark  D.  Fisk 

Mission  Research  Corporation 

735  State  Street 

P.O  Drawer  7 19 

Santa  Barbara,  CA  93102 

Prof  Stanley  Flatte 
Applied  Sciences  Building 
l  mveisity  of  California,  Santa  Cruz 
Santa  Cruz.  CA  95064 


Dr.  John  Foley 
NEK-Geo  Sciences 
1 1  (X)  Crown  Colony  Drive 
Quincy.  MA  02169 


Prof.  Donald  Forsyth 
Department  of  Geological  Sciences 
Brown  University 
Providence.  RI  02912 


Dr.  Art  Frunkel 
U.S.  Geological  Survey 
922  National  Center 
Reston,  VA  22092 


Dr.  Cliff  Frolich 
Institute  of  Geophysics 
8701  North  Mopac 
Austin,  TX  78759 


Dr.  Holly  Given 
1GPP,  A-025 

Scripps  Institute  of  Oceanography 
University  of  California.  San  Diego 
La  Jolla,  CA  92093 

Dr.  Jeffrev  W.  Given 
SAiC 

10260  Campus  Point  Drive 
San  Diego,  CA  92121 


Dr.  Dale  Glover 
Defense  Intelligence  Attencv 
A'lTN:  ODT-1B 
Washington,  DC  20301 


Dr.  Indra  Gupta 
Teledyne  Geotech 
314  Montgomery  Street 
Alexanderia,  VA  22314 


Dan  N.  Hagedon 
Pacific  Northwest  Laboratories 
Battelle  Boulevard 
Richland,  WA  99352 


Dr.  James  Hannon 

Lawrence  Livermore  National  Laboratory 

P.O.  Box  808 

L-205 

Livermore,  CA  94550 

Dr.  Roger  Hansen 
HQ  APTAC/TTR 
1 30  South  Highway  A 1 A 
Patrick  AFB,  FL  32925-3002 


Prof.  David  G.  Harkrider 
Seismological  Laboratory 
Division  of  Geological  &  Planetary  Scieno 
California  Institute  of  Technology 
Pasadena,  CA  91125 

Prof.  Danny  Harvcv 
CIRES 

University  of  Colorado 
Boulder,  CO  80309 


Prof.  Donald  V.  Helmberger 
Seismological  Laboratory 
Division  of  Geological  &  Planetary  Sciences 
California  Institute  of  Technology 
Pasadena,  CA  91125 

Prof.  Eugene  Herrin 

Institute  for  the  Study  of  Earth  and  Man 

Geophysical  Laboratory 

Southern  Methodist  University 

Dallas,  TX  75275 

Prof.  Robert  B.  Herrmann 

Department  of  Earth  &  Atmospheric  Sciences 

St.  Louis  University 

St.  Louis,  MO  63156 


Prof.  Lane  R.  Johnson 
Seismographic  Station 
University  of  California 
Berkeley,  CA  94720 


Prof.  Thomas  H.  Jordan 
Department  of  Earth,  Atmospheric  & 
Planetary  Sciences 

Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 

Prof.  Alan  Kafka 

Department  of  Geology  &  Geophysics 
Boston  College 
Chestnut  Hill,  MA  02167 


Robert  C.  Kemerait 
ENSCO,  Inc. 

445  Pineda  Court 
Melbourne,  FL  32940 


Dr.  Karl  Koch 

Institute  for  the  Study  of  Earth  and  Man 
Geophysical  Laboratory 
Southern  Methodist  University 
Dallas,  Tx  75275 

Dr.  Max  Koontz 
U.S.  Dept,  of  Energy/DP  5 
Forres tal  Building 
1000  Independence  Avenue 
Washington,  DC  20585 

Dr.  Richard  LaCoss 

MIT  Lincoln  Laboratory,  M-200B 

P.O.  Box  73 

Lexington,  MA  02173-0073 


Dr.  Fred  K  Lamb 

University  of  Illinois  at  Urbana-C'hampaign 
Department  of  Physics 
1110  West  Green  Street 
Urbana.  II.  61, SOI 

Prof.  Charles  A.  Langston 
Geosciences  Department 
403  Deike  Building 
The  Pennsylvania  State  University 
University  Park.  PA  16802 

Jim  Law  son,  Chief  Geophy  sicist 
Oklahoma  Geological  Survey 
Oklahoma  Geophysical  Observatory 
P.O.  Box  S 

Leonard.  OK  74043-0008 

Prof.  Thorne  Lay 

Institute  of  Tectonics 

Earth  Science  Board 

University  of  California,  Santa  Cru/ 

Santa  Cruz,  CA  95064 

Dr.  Wilham  l^eith 
U.S.  Geological  Survey 
Mail  Stop  928 
Reston,  VA  22092 


Mr  James  F.  Lewkowicz 
Phillips  Laboratory/GPEH 
29  Randolph  Road 

Hunscom  AFB.  MA  01731-30KH  2  copies 


Mr.  Alfred  Lieberman 

ACDA/VI-OA  State  Department  Building 

Room  5726 

320-2 1st  Street,  NW 

Washington,  DC  20451 

Prof.  L.  Timothy  Long 
School  of  Geophysical  Sciences 
Georgia  Institute  of  Technology 
Atlanta,  GA  30332 


Dr.  Randolph  Martin.  Ill 
New  England  Research,  Inc. 

76  Olcott  Drive 

White  River  Junction,  VT  05(X)1 


Dr,  Robert  Masse 
Denver  Federal  Building 
Box  25046.  Mail  Stop  907 
Denver.  CO  80225 


3 


Dr,  Gary  McCarter 
Department  of  Physics 
Southern  Methodist  University 
Dallas.  TX  75275 


Prof.  Thomas  V.  McEvilly 
Seismographic  Station 
Universitv  of  California 
Berkeley,' CA  94720 


Dr.  An  McGarr 
U  S.  Geological  Survey 
Mail  Stop  977 
L’.S.  Geological  Survey 
Menlo  Park.  CA  94025 

Dr.  Keith  L.  McUaughlin 
S-CUBED 

A  Division  of  Maxwell  Laboratory 

P.O.  Box  1620 

La  Jolla.  CA  92038-1620 

Stephen  Miller  &  Dr.  Alexander  Florence 

SRI  International 

333  Ravenswood  Avenue 

Box  AF  1  16 

Menlo  Park.  CA  94025-3493 

Prof.  Bernard  Minster 
1GPP,  A-025 

Scripps  Institute  of  Oceanography 
Universitv  of  California,  San  Dieeo 
La  Jolla.  CA  92093 

Prof.  Brian  J.  Mitchell 

Department  of  Earth  &  Atmospheric  Sciences 
St  1  otiis  University 
St.  Louis,  MO  63156 


Mr.  Jack  Murphv 
S-CUBED 

A  Division  of  Maxwell  Laboratory 
1 1 800  Sunrise  Valley  Drive,  Suite  1212 
Heston.  V A  22091  (2  Copies) 

Dr.  Keith  K.  Nakanishi 

Lawrence  Livermore  National  Laboratory 

L-025 

P.O.  Box  80S 
Livermore,  CA  94550 

Dr.  Carl  Newton 

Los  Alamo.>  National  I  'horutory 

P.O.  Box  1663 

Mail  Stop  (035,  Group  ESS-3 
Los  Alamos.  NM  87545 


Dr.  Bao  Nguven 
HQ  AfTAC/TTR 
130  South  Highway  A1 A 
Patrick  AFB,  FL  32925-3002 


Prof.  John  A.  Orctitt 
IGPP.  A-025 

Scripps  Institute  of  (Oceanography 
University  of  California.  San  Diego 
La  Jolla,  CA  92093 

Prof.  Jeffrey  Park 

Kline  Geology  Laboratory 

P.O.  Box  6666 

New  Haven,  CT  0(.M  1-8130 


Dr.  Howard  Patton 

Lawrence  Livermore  National  Laboratory 
L-025 

P  O.  Box  808 
Livermore,  CA  94550 

Dr.  Frank  Pi  lotto 
HQ  AFTAC/TT 
130  South  Highway  A 1 A 
Patrick  AFB,  FL  32925-3002 


Dr.  Jay  J.  Puili 
Radix  Systems,  Inc. 

201  Perry  Parkway 
Gaithersburg,  MD  20877 


Dr.  Robert  Reinke 
A  TUN:  FCTVTD 
Field  Command 
Defense  Nuclear  Agency 
Kirtland  AFB,  NM"  87il5 

Prof.  Paul  G.  Richards 
Lamont-Doherty  Geological  Observatory 
of  Columbia  Universitv 
Palisades.  NY  10964 


Mr.  Wilmer  Rivers 
Teledyne  Geotech 
314  Montgomery  Street 
Alexandria,  VA  22314 


Dr.  George  Rothe 
HQ  AFTAQTTR 
130  South  Highway  A 1 A 
Patrit  k  AFB,  hi.  32925  3002 


4 


Dr.  Alan  S.  Rvall,  Jr. 
DARPA/NMRO 
3701  North  Fairfax  Drive 
Arlington.  VA  22209-1714 


Dr.  Richard  Sailor 
TASC,  Inc. 

55  Walkers  Brook  Drive 
Reading,  MA  01867 


Prof.  Charles  G.  Sammis 
Center  for  Earth  Sciences 
University  of  Southern  California 
University  Park 
Los  Angeles,  CA  90089-0741 

Prof.  Christopher  H.  Scholz 
Lamont-Doherty  Geological  Observatory 
of  Columbia  University 
Palisades.  NY  10964 


Dr.  Susan  Schwartz 
Institute  of  Tectonics 
1 156  High  Street 
Santa  Cruz,  CA  95064 


Secretary  of  the  Air  Force 
(SAFRD) 

Washington,  DC  20330 


Office  of  the  Secretary  of  Defense 
DDR&E 

Washington,  DC  20330 


Thomas  J.  Sereno,  Jr. 

Science  Application  Int'l  Corp. 
10260  Campus  Point  Drive 
San  Diego,  CA  92121 


Dr.  Michael  Shore 
Defense  Nuclear  Agency/SPSS 
6801  Telegraph  Road 
Alexandria,  VA  22310 


Dr.  Robert  Shumway 
University  of  California  Davis 
Division  of  Statistics 
Davis,  CA  95616 


Dr.  Matthew  Si  hoi 
Virginia  Tech 

Seismological  Observatory 
4044  I  Vrrine  I  tail 
Blacksburg.' V A  24001 0420 

Prof.  David  G.  Simpson 
IRIS.  Inc. 

1616  North  Fort  Mver  Drive 
Suite  1050 
Arlington.  VA  22209 

Donald  L.  Springer 

Lawrence  Livermore  National  Lahor.uor. 

L  025 

P.O.  Box  SOS 
Livermore.  ('A  94550 

Dr.  Jeff  rev  Stevens 
S-CUBFD 

A  Division  of  Maxwell  Laboratorv 
P.O.  Box  1620 
La  Jolla.  CA  92038-1620 

Lt.  Col.  Jim  Stobie 
ATTN:  AFOSR/NL 
1 10  Duncan  Avenue 
Bolling  AFB 

Washington,  DC  20332-0001 
Prof.  Brian  Stump 

Institute  for  the  Study  of  Earth  <N  Man 
Geophysical  Laboratory 
Southern  Methodist  University 
Dallas.  TX  75275 

Prof.  Jeremiah  Sullivan 

University  of  Illinois  at  Urbana-Chamnaign 

Department  of  Physics 

1110  West  Green  Street 

Urbana,  1L  61801 

Prof.  L.  Sykes 

Lamont-Doherty  Geological  Observatory 
of  Columbia  University 
Palisades,  NY  10964 


Dr.  David  Tavlor 
ENSCO.  Inc' 

445  Pineda  Court 
Melbourne,  FL  32940 


Dr.  Steven  R.  Taylor 

Los  Alamos  National  Laboratorv 

P.O.  Box  1663 

Mail  Stop  C335 

Los  Alamos,  NM  87545 


5 


Prof.  Clifford  Thurber 
University  of  Wisconsin-Madison 
Department  of  Geology  &  Geophysics 
1215  West  Dayton  Street 
Madison,  WS  53706 

Prof.  M.  Nafi  Toksoz 
Earth  Resources  Lab 
Massachusetts  Institute  of  Technology 
42  Carleton  Street 
Cambridge,  MA  02142 

Dr  Larry  Turnbull 
CIA-OSWR/NED 
Washington,  DC  20505 


DARPA/PM 

3701  North  Fairfax  Drive 
Arlington,  V  A  22203-1714 


DARPA/RMO/kETRlEVAL 
3701  North  Fairfax  Drive 
Arlington,  V A  22203-1714 


DARPA/RMO/SECURITY  OFFICE 
3701  North  Fairfax  Drive 
Arlington,  V A  22203-1714 


Dr.  Gregory  van  der  Vink 
IRIS,  Inc. 

1616  North  Fort  Myer  Drive 
Suite  1050 

Arlington,  VA  22209 

Dr.  Karl  Veith 
EG&G 

5211  Auth  Road 
Suite  240 

Suitland,  MD  20746 

Prof.  Terry  C.  Wallace 
Department  of  Geosciences 
Building  #77 
University  of  Arizona 
Tuscon,  AZ  85721 

Dr.  Thomas  Weaver 

I  . os  Alamos  National  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335 

Los  Alamos,  NM  87545 


HQ  DNA 

ATTN:  Technical  Library 
Washington,  DC  20305 


Defense  Intelligence  Agency 

Directorate  for  Scientific  &  Technical  Intelligence 

ATTN:  DTIB 

Washington,  DC  20340-6158 


Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  V A  22314  (2  Copies) 


TACTEC 

Battelle  Memorial  Institute 
505  King  Avenue 

Columbus,  OH  43201  (Final  Report) 


Dr.  William  Wortman 
Mission  Research  Corporation 
8560  Cinderbed  Road 
Suite  700 

Newington,  VA  22122 

Prof.  Francis  T.  Wu 
Department  of  Geological  Sciences 
State  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 

AFrAC/CA 

(STINFO) 

Patrick  AFB,  FL  32925-6001 


Phillips  Laboratory 

ATTN.  XPG 

29  Randolph  Road 

Hanscom  AFB,  MA  01731  -3010 


Phillips  Laboratory 

ATTN:  GPE 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 


Phillips  Laboratory 
ATTN:  TSML 
5  Wright  Street 

Hanscom  AFB,  MA  01731-3004 


6 


Phillips  Laboratory 

ATTN:  PL/SUL 

3550  Aberdeen  Ave  SE 

Kirtlarul,  NM  871 17-5776  (2  copies) 


Dr.  Svein  Mykkeltveit 
NTNT/NORSAR 
P.O.  Box  5  I 

N  2007  Kjcllcr,  NORWAY  (3  Copies j 


Dr.  Michel  Bouchon 

I.R.I.G.M.-B.P.  68 

38402  St.  Martin  D’Heres 

Cedex,  FRANCE 

Prof.  Keith  Priestley 

University  of  Cambridge 

Bullard  Labs,  Dept,  of  Earth  Sciences 
Madinglev  Rise,  Madingley  Road 

Cambridge  CB3  OEZ,  ENGLAND 

Dr.  Michel  Campillo 

Observatoire  de  Grenoble 

I.R.I.G.M.-B.P.  53 

38041  Grenoble,  FRANCE 

Dr.  Jorg  Schiittenhardt 

Federal  Institute  for  Geosciences  &  Nat  !  Re 
Postfach  510153 

D-3000  Hannover  5 1,  GERMAN Y 

Dr.  Kin  Yip  Chun 

Geophysics  Division 

Physics  Department 

University  of  Toronto 

Ontario,  CANADA 

Dr.  Johannes  Schweitzer 

Institute  of  Geophysics 

Ruhr  University/Bochum 

P.O.  Box  1  102148 

4360  Bochum  1 ,  GERMANY 

Prof.  Hans- Peter  Harjes 

Institute  for  Geophysic 

Ruhr  University/Bochum 

P.O.  Box  102148 

4630  Bochum  1,  GERMANY 

Trust  &  Verifv 

VERTIC 

8  John  Adam  Street 

London  WC2N  6EZ,  ENGLAND 

Prof.  Eystein  Husebye 

NTNF/NORSAR 

P.O.  Box  51 

N-2007  Kjeller,  NORWAY 

David  Jepsen 

Acting  Head,  Nuclear  Monitoring  Section 

Bureau  of  Mineral  Resources 

Geology  and  Geophysics 

G.P.O.  Box  378,  Canberra,  AUSTRALIA 

Ms.  Eva  Johannisson 
Senior  Research  Officer 
FOA 

S-172  90  Sundbyberg,  SWEDEN 


Dr.  Peter  Marshall 
Procurement  Executive 
Ministry  of  Defense 
Blacknest,  Brimpton 

Reading  FG7-FRS,  UNITED  KINGDOM 

Dr.  Bernard  Massinon,  Dr.  Pierre  Mechler 
Societe  Radiomana 
27  rue  Claude  Bernard 
75005  Paris,  FRANCE  (2  Copies) 


7 


